diff --git a/THIRD_PARTY_LICENSES/MT19937-64-LICENSE b/THIRD_PARTY_LICENSES/MT19937-64-LICENSE new file mode 100644 index 00000000..2642951c --- /dev/null +++ b/THIRD_PARTY_LICENSES/MT19937-64-LICENSE @@ -0,0 +1,52 @@ +A C-program for MT19937-64 (2004/9/29 version). +Coded by Takuji Nishimura and Makoto Matsumoto. + +This is a 64-bit version of Mersenne Twister pseudorandom number +generator. + +Before using, initialize the state by using init_genrand64(seed) +or init_by_array64(init_key, key_length). + +Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura, +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +References: +T. Nishimura, ``Tables of 64-bit Mersenne Twisters'' + ACM Transactions on Modeling and + Computer Simulation 10. (2000) 348--357. +M. Matsumoto and T. Nishimura, + ``Mersenne Twister: a 623-dimensionally equidistributed + uniform pseudorandom number generator'' + ACM Transactions on Modeling and + Computer Simulation 8. (Jan. 1998) 3--30. + +Any feedback is very welcome. +http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html +email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces) diff --git a/THIRD_PARTY_LICENSES/MT19937-LICENSE b/THIRD_PARTY_LICENSES/MT19937-LICENSE new file mode 100644 index 00000000..08207eb5 --- /dev/null +++ b/THIRD_PARTY_LICENSES/MT19937-LICENSE @@ -0,0 +1,40 @@ +A C-program for MT19937, with initialization improved 2002/1/26. +Coded by Takuji Nishimura and Makoto Matsumoto. + +Before using, initialize the state by using init_genrand(seed) +or init_by_array(init_key, key_length). + +Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + +Any feedback is very welcome. +http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html +email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) diff --git a/mathext/prng/mt19937.go b/mathext/prng/mt19937.go new file mode 100644 index 00000000..bf375bf8 --- /dev/null +++ b/mathext/prng/mt19937.go @@ -0,0 +1,151 @@ +// Copyright ©2019 The Gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +// Original C program copyright Takuji Nishimura and Makoto Matsumoto 2002. +// http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c + +package prng + +import ( + "encoding/binary" + "io" +) + +const ( + mt19937N = 624 + mt19937M = 397 + mt19937matrixA = 0x9908b0df + mt19937UpperMask = 0x80000000 + mt19937LowerMask = 0x7fffffff +) + +// MT19937 implements the 32 bit Mersenne Twister PRNG. MT19937 +// is the default PRNG for a wide variety of programming systems. +// See https://en.wikipedia.org/wiki/Mersenne_Twister. +type MT19937 struct { + mt [mt19937N]uint32 + mti uint32 +} + +// NewMT19937 returns a new MT19937 PRNG. The returned PRNG will +// use the default seed 5489 unless the Seed method is called with +// another value. +func NewMT19937() *MT19937 { + return &MT19937{mti: mt19937N + 1} +} + +// Seed uses the provided seed value to initialize the generator to a +// deterministic state. Only the lower 32 bits of seed are used to seed +// the PRNG. +func (src *MT19937) Seed(seed uint64) { + src.mt[0] = uint32(seed) + for src.mti = 1; src.mti < mt19937N; src.mti++ { + src.mt[src.mti] = (1812433253*(src.mt[src.mti-1]^(src.mt[src.mti-1]>>30)) + src.mti) + } +} + +// SeedFromKeys uses the provided seed key value to initialize the +// generator to a deterministic state. It is provided for compatibility +// with C implementations. +func (src *MT19937) SeedFromKeys(keys []uint32) { + src.Seed(19650218) + i := uint32(1) + j := uint32(0) + k := uint32(mt19937N) + if k <= uint32(len(keys)) { + k = uint32(len(keys)) + } + for ; k != 0; k-- { + src.mt[i] = (src.mt[i] ^ ((src.mt[i-1] ^ (src.mt[i-1] >> 30)) * 1664525)) + keys[j] + j // Non linear. + i++ + j++ + if i >= mt19937N { + src.mt[0] = src.mt[mt19937N-1] + i = 1 + } + if j >= uint32(len(keys)) { + j = 0 + } + } + for k = mt19937N - 1; k != 0; k-- { + src.mt[i] = (src.mt[i] ^ ((src.mt[i-1] ^ (src.mt[i-1] >> 30)) * 1566083941)) - i // Non linear. + i++ + if i >= mt19937N { + src.mt[0] = src.mt[mt19937N-1] + i = 1 + } + } + src.mt[0] = 0x80000000 // MSB is 1; assuring non-zero initial array. +} + +// Uint32 returns a pseudo-random 32-bit unsigned integer as a uint32. +func (src *MT19937) Uint32() uint32 { + mag01 := [2]uint32{0x0, mt19937matrixA} + + var y uint32 + if src.mti >= mt19937N { // Generate mt19937N words at one time. + if src.mti == mt19937N+1 { + // If Seed() has not been called + // a default initial seed is used. + src.Seed(5489) + } + + var kk int + for ; kk < mt19937N-mt19937M; kk++ { + y = (src.mt[kk] & mt19937UpperMask) | (src.mt[kk+1] & mt19937LowerMask) + src.mt[kk] = src.mt[kk+mt19937M] ^ (y >> 1) ^ mag01[y&0x1] + } + for ; kk < mt19937N-1; kk++ { + y = (src.mt[kk] & mt19937UpperMask) | (src.mt[kk+1] & mt19937LowerMask) + src.mt[kk] = src.mt[kk+(mt19937M-mt19937N)] ^ (y >> 1) ^ mag01[y&0x1] + } + y = (src.mt[mt19937N-1] & mt19937UpperMask) | (src.mt[0] & mt19937LowerMask) + src.mt[mt19937N-1] = src.mt[mt19937M-1] ^ (y >> 1) ^ mag01[y&0x1] + + src.mti = 0 + } + + y = src.mt[src.mti] + src.mti++ + + // Tempering. + y ^= (y >> 11) + y ^= (y << 7) & 0x9d2c5680 + y ^= (y << 15) & 0xefc60000 + y ^= (y >> 18) + + return y +} + +// Uint64 returns a pseudo-random 64-bit unsigned integer as a uint64. +// It makes use of two calls to Uint32 placing the first result in the +// upper bits and the second result in the lower bits of the returned +// value. +func (src *MT19937) Uint64() uint64 { + h := uint64(src.Uint32()) + l := uint64(src.Uint32()) + return h<<32 | l +} + +// MarshalBinary returns the binary representation of the current state of the generator. +func (src *MT19937) MarshalBinary() ([]byte, error) { + var buf [(mt19937N + 1) * 4]byte + for i := 0; i < mt19937N; i++ { + binary.BigEndian.PutUint32(buf[i*4:(i+1)*4], src.mt[i]) + } + binary.BigEndian.PutUint32(buf[mt19937N*4:], src.mti) + return buf[:], nil +} + +// UnmarshalBinary sets the state of the generator to the state represented in data. +func (src *MT19937) UnmarshalBinary(data []byte) error { + if len(data) < (mt19937N+1)*4 { + return io.ErrUnexpectedEOF + } + for i := 0; i < mt19937N; i++ { + src.mt[i] = binary.BigEndian.Uint32(data[i*4 : (i+1)*4]) + } + src.mti = binary.BigEndian.Uint32(data[mt19937N*4:]) + return nil +} diff --git a/mathext/prng/mt19937_64.go b/mathext/prng/mt19937_64.go new file mode 100644 index 00000000..19592211 --- /dev/null +++ b/mathext/prng/mt19937_64.go @@ -0,0 +1,142 @@ +// Copyright ©2019 The Gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +// Original C program copyright Takuji Nishimura and Makoto Matsumoto 2004. +// http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/mt19937-64.c + +package prng + +import ( + "encoding/binary" + "io" +) + +const ( + mt19937_64NN = 312 + mt19937_64MM = 156 + mt19937_64MatrixA = 0xB5026F5AA96619E9 + mt19937_64UpperMask = 0xFFFFFFFF80000000 + mt19937_64LowerMask = 0x7FFFFFFF +) + +// MT19937_64 implements the 64 bit Mersenne Twister PRNG. MT19937_64 +// is the 64 bit version of MT19937, it has the same sized state, but +// generates a different sequence. +// See https://en.wikipedia.org/wiki/Mersenne_Twister. +type MT19937_64 struct { + mt [mt19937_64NN]uint64 + mti uint64 +} + +// NewMT19937_64 returns a new MT19937_64 PRNG. The returned PRNG will +// use the default seed 5489 unless the Seed method is called with +// another value. +func NewMT19937_64() *MT19937_64 { + return &MT19937_64{mti: mt19937_64NN + 1} +} + +// Seed uses the provided seed value to initialize the generator to a +// deterministic state. +func (src *MT19937_64) Seed(seed uint64) { + src.mt[0] = seed + for src.mti = 1; src.mti < mt19937_64NN; src.mti++ { + src.mt[src.mti] = (6364136223846793005*(src.mt[src.mti-1]^(src.mt[src.mti-1]>>62)) + src.mti) + } +} + +// SeedFromKeys uses the provided seed key value to initialize the +// generator to a deterministic state. It is provided for compatibility +// with C implementations. +func (src *MT19937_64) SeedFromKeys(keys []uint64) { + src.Seed(19650218) + i := uint64(1) + j := uint64(0) + k := uint64(mt19937_64NN) + if k <= uint64(len(keys)) { + k = uint64(len(keys)) + } + for ; k != 0; k-- { + src.mt[i] = (src.mt[i] ^ ((src.mt[i-1] ^ (src.mt[i-1] >> 62)) * 3935559000370003845)) + keys[j] + j // Non linear. + i++ + j++ + if i >= mt19937_64NN { + src.mt[0] = src.mt[mt19937_64NN-1] + i = 1 + } + if j >= uint64(len(keys)) { + j = 0 + } + } + for k = mt19937_64NN - 1; k != 0; k-- { + src.mt[i] = (src.mt[i] ^ ((src.mt[i-1] ^ (src.mt[i-1] >> 62)) * 2862933555777941757)) - i // Non linear. + i++ + if i >= mt19937_64NN { + src.mt[0] = src.mt[mt19937_64NN-1] + i = 1 + } + } + + src.mt[0] = 1 << 63 /* MSB is 1; assuring non-zero initial array */ +} + +// Uint64 returns a pseudo-random 64-bit unsigned integer as a uint64. +func (src *MT19937_64) Uint64() uint64 { + mag01 := [2]uint64{0, mt19937_64MatrixA} + + var x uint64 + if src.mti >= mt19937_64NN { // Generate mt19937_64NN words at one time. + if src.mti == mt19937_64NN+1 { + // If Seed() has not been called + // a default initial seed is used. + src.Seed(5489) + } + + var i int + for ; i < mt19937_64NN-mt19937_64MM; i++ { + x = (src.mt[i] & mt19937_64UpperMask) | (src.mt[i+1] & mt19937_64LowerMask) + src.mt[i] = src.mt[i+mt19937_64MM] ^ (x >> 1) ^ mag01[(int)(x&0x1)] + } + for ; i < mt19937_64NN-1; i++ { + x = (src.mt[i] & mt19937_64UpperMask) | (src.mt[i+1] & mt19937_64LowerMask) + src.mt[i] = src.mt[i+(mt19937_64MM-mt19937_64NN)] ^ (x >> 1) ^ mag01[(int)(x&0x1)] + } + x = (src.mt[mt19937_64NN-1] & mt19937_64UpperMask) | (src.mt[0] & mt19937_64LowerMask) + src.mt[mt19937_64NN-1] = src.mt[mt19937_64MM-1] ^ (x >> 1) ^ mag01[(int)(x&0x1)] + + src.mti = 0 + } + + x = src.mt[src.mti] + src.mti++ + + // Tempering. + x ^= (x >> 29) & 0x5555555555555555 + x ^= (x << 17) & 0x71D67FFFEDA60000 + x ^= (x << 37) & 0xFFF7EEE000000000 + x ^= (x >> 43) + + return x +} + +// MarshalBinary returns the binary representation of the current state of the generator. +func (src *MT19937_64) MarshalBinary() ([]byte, error) { + var buf [(mt19937_64NN + 1) * 8]byte + for i := 0; i < mt19937_64NN; i++ { + binary.BigEndian.PutUint64(buf[i*8:(i+1)*8], src.mt[i]) + } + binary.BigEndian.PutUint64(buf[mt19937_64NN*8:], src.mti) + return buf[:], nil +} + +// UnmarshalBinary sets the state of the generator to the state represented in data. +func (src *MT19937_64) UnmarshalBinary(data []byte) error { + if len(data) < (mt19937_64NN+1)*8 { + return io.ErrUnexpectedEOF + } + for i := 0; i < mt19937_64NN; i++ { + src.mt[i] = binary.BigEndian.Uint64(data[i*8 : (i+1)*8]) + } + src.mti = binary.BigEndian.Uint64(data[mt19937_64NN*8:]) + return nil +} diff --git a/mathext/prng/mt19937_64_test.go b/mathext/prng/mt19937_64_test.go new file mode 100644 index 00000000..3e1c8222 --- /dev/null +++ b/mathext/prng/mt19937_64_test.go @@ -0,0 +1,88 @@ +// Copyright ©2019 The Gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package prng + +import ( + "testing" + "time" + + "golang.org/x/exp/rand" +) + +var _ rand.Source = (*MT19937_64)(nil) + +// Random values in tests are produced by 40 iterations of the C code +// with or without an initial seed array. + +func TestMT19937_64(t *testing.T) { + want := []uint64{ + 14514284786278117030, 4620546740167642908, 13109570281517897720, 17462938647148434322, 355488278567739596, + 7469126240319926998, 4635995468481642529, 418970542659199878, 9604170989252516556, 6358044926049913402, + 5058016125798318033, 10349215569089701407, 2583272014892537200, 10032373690199166667, 9627645531742285868, + 15810285301089087632, 9219209713614924562, 7736011505917826031, 13729552270962724157, 4596340717661012313, + 4413874586873285858, 5904155143473820934, 16795776195466785825, 3040631852046752166, 4529279813148173111, + 3658352497551999605, 13205889818278417278, 17853215078830450730, 14193508720503142180, 1488787817663097441, + 8484116316263611556, 4745643133208116498, 14333959900198994173, 10770733876927207790, 17529942701849009476, + 8081518017574486547, 5945178879512507902, 9821139136195250096, 4728986788662773602, 840062144447779464, + } + + mt := NewMT19937_64() + for i := range want { + got := mt.Uint64() + if got != want[i] { + t.Errorf("unexpected random value at iteration %d: got:%d want:%d", i, got, want[i]) + } + } +} + +func TestMT19937_64SeedFromKeys(t *testing.T) { + want := []uint64{ + 7266447313870364031, 4946485549665804864, 16945909448695747420, 16394063075524226720, 4873882236456199058, + 14877448043947020171, 6740343660852211943, 13857871200353263164, 5249110015610582907, 10205081126064480383, + 1235879089597390050, 17320312680810499042, 16489141110565194782, 8942268601720066061, 13520575722002588570, + 14226945236717732373, 9383926873555417063, 15690281668532552105, 11510704754157191257, 15864264574919463609, + 6489677788245343319, 5112602299894754389, 10828930062652518694, 15942305434158995996, 15445717675088218264, + 4764500002345775851, 14673753115101942098, 236502320419669032, 13670483975188204088, 14931360615268175698, + 8904234204977263924, 12836915408046564963, 12120302420213647524, 15755110976537356441, 5405758943702519480, + 10951858968426898805, 17251681303478610375, 4144140664012008120, 18286145806977825275, 13075804672185204371, + } + + mt := NewMT19937_64() + mt.SeedFromKeys([]uint64{0x12345, 0x23456, 0x34567, 0x45678}) + for i := range want { + got := mt.Uint64() + if got != want[i] { + t.Errorf("unexpected random value at iteration %d: got:%d want:%d", i, got, want[i]) + } + } +} + +func TestMT19937_64RoundTrip(t *testing.T) { + var src MT19937_64 + src.Seed(uint64(time.Now().Unix())) + + src.Uint64() // Step PRNG once to makes sure states are mixed. + + buf, err := src.MarshalBinary() + if err != nil { + t.Errorf("unexpected error marshaling state: %v", err) + } + + var dst MT19937_64 + // Get dst into a non-zero state. + dst.Seed(1) + for i := 0; i < 10; i++ { + dst.Uint64() + } + + err = dst.UnmarshalBinary(buf) + if err != nil { + t.Errorf("unexpected error unmarshaling state: %v", err) + } + + if dst != src { + t.Errorf("mismatch between generator states: got:%+v want:%+v", dst, src) + } +} diff --git a/mathext/prng/mt19937_test.go b/mathext/prng/mt19937_test.go new file mode 100644 index 00000000..e63e67fe --- /dev/null +++ b/mathext/prng/mt19937_test.go @@ -0,0 +1,88 @@ +// Copyright ©2019 The Gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package prng + +import ( + "testing" + "time" + + "golang.org/x/exp/rand" +) + +var _ rand.Source = (*MT19937)(nil) + +// Random values in tests are produced by 40 iterations of the C code +// with or without an initial seed array. + +func TestMT19937(t *testing.T) { + want := []uint32{ + 3499211612, 581869302, 3890346734, 3586334585, 545404204, + 4161255391, 3922919429, 949333985, 2715962298, 1323567403, + 418932835, 2350294565, 1196140740, 809094426, 2348838239, + 4264392720, 4112460519, 4279768804, 4144164697, 4156218106, + 676943009, 3117454609, 4168664243, 4213834039, 4111000746, + 471852626, 2084672536, 3427838553, 3437178460, 1275731771, + 609397212, 20544909, 1811450929, 483031418, 3933054126, + 2747762695, 3402504553, 3772830893, 4120988587, 2163214728, + } + + mt := NewMT19937() + for i := range want { + got := mt.Uint32() + if got != want[i] { + t.Errorf("unexpected random value at iteration %d: got:%d want:%d", i, got, want[i]) + } + } +} + +func TestMT19937SeedFromKeys(t *testing.T) { + want := []uint32{ + 1067595299, 955945823, 477289528, 4107218783, 4228976476, + 3344332714, 3355579695, 227628506, 810200273, 2591290167, + 2560260675, 3242736208, 646746669, 1479517882, 4245472273, + 1143372638, 3863670494, 3221021970, 1773610557, 1138697238, + 1421897700, 1269916527, 2859934041, 1764463362, 3874892047, + 3965319921, 72549643, 2383988930, 2600218693, 3237492380, + 2792901476, 725331109, 605841842, 271258942, 715137098, + 3297999536, 1322965544, 4229579109, 1395091102, 3735697720, + } + + mt := NewMT19937() + mt.SeedFromKeys([]uint32{0x123, 0x234, 0x345, 0x456}) + for i := range want { + got := mt.Uint32() + if got != want[i] { + t.Errorf("unexpected random value at iteration %d: got:%d want:%d", i, got, want[i]) + } + } +} + +func TestMT19937RoundTrip(t *testing.T) { + var src MT19937 + src.Seed(uint64(time.Now().Unix())) + + src.Uint64() // Step PRNG once to makes sure states are mixed. + + buf, err := src.MarshalBinary() + if err != nil { + t.Errorf("unexpected error marshaling state: %v", err) + } + + var dst MT19937 + // Get dst into a non-zero state. + dst.Seed(1) + for i := 0; i < 10; i++ { + dst.Uint64() + } + + err = dst.UnmarshalBinary(buf) + if err != nil { + t.Errorf("unexpected error unmarshaling state: %v", err) + } + + if dst != src { + t.Errorf("mismatch between generator states: got:%+v want:%+v", dst, src) + } +}