From b631207ca737a6b1e3bd30d10d687ccaf9b4918b Mon Sep 17 00:00:00 2001 From: Kirk Haines Date: Thu, 5 Dec 2024 10:05:43 +0100 Subject: [PATCH] feat(examples): Add a useful set of high quality pseudo-random number generators (#2868) I ported a number of my pseudo-random number generator implementations from Ruby to gno while traveling to the retreat last weekend as an exercise in expanding my comfort level with gno code, and expanding my understanding of some of the code internals, while contributing code that others may find interesting or useful. I added two xorshift generators, xorshift64* and xorshiftr128+. These are both many times faster than the PCG generator that is the gno default, and produce high quality randomness with great statistical qualities. In addition to these, I added both the 32-bit ISAAC implementation (with an added function to return 64 bit values), and the 64-bit ISAAC implementation. ISAAC is a stellar pseudo-random number generator. Both implementations are significantly faster than PCG (though not near so fast as the xorshift algorithms), while producing extremely high quality, cryptographically secure randomness that can not be differentiated from real randomness. All of these were built to be compatible with the standard Rand() implementation. This means that any of these can be used as a drop-in replacement for the default PCG algorithm: ``` source = isaac.New() prng := rand.New(source) ``` All of these leverage the `gno.land/p/demo/entropy` package to assist with seeding if no seed is provided. In the case of the ISAAC algorithms, they require 256 uint values for their seed, so they leverage a combination of `entropy` and `xorshiftr128+` to generate any missing numbers in the provided seed. I also added a function to entropy to return uint64, to facilitate using it for seeding. I added tests to entropy, and wrote tests for the other generators, as well. There are a few other things that ended up in this PR. In order to make some fact based assertions about the performance of these generators, I included some code that can be ran via `gno run -expr`. i.e. `gno run -expr 'averageISAAC()' isaac.gno` that can be used to get some benchmarks and some very simple self-statistical-analysis on the results, and when I did so, I discovered that the current `ufmt.Sprintf` implementation didn't support any of the float output flags. I added float support to it's capabilities, which, in turn, required adding `FormatFloat` to the `strconv.gno/strconv.go` implementation in the standard library. I added a test to cover this. I also noticed that there is a test in `tm2/pkg/p2p` that is failing on both master and my branch. Specifically, there is a call to `sw.Logger.Error()` that passes a message and an error, but not `"err"` before the error. Adding that seemed to clear up the build failure. This, specifically, is line 222 of `switch.go`. Currently there is one failing test, which is the code coverage check on tm2, because it is non-obvious to me how to setup a test to properly exercise that one changed line.
Contributors' checklist... - [X] Added new tests, or not needed, or not feasible - [X] Provided an example (e.g. screenshot) to aid review or the PR is self-explanatory - [X] Updated the official documentation or not needed - [X] No breaking changes were made, or a `BREAKING CHANGE: xxx` message was included in the description - [X] Added references to related issues and PRs
--------- Co-authored-by: Morgan Co-authored-by: Nathan Toups <612924+n2p5@users.noreply.github.com> Co-authored-by: Morgan --- examples/gno.land/p/demo/entropy/entropy.gno | 8 + .../gno.land/p/demo/entropy/entropy_test.gno | 32 ++ .../gno.land/p/demo/entropy/z_filetest.gno | 6 + examples/gno.land/p/demo/ufmt/ufmt_test.gno | 6 + .../gno.land/p/wyhaines/rand/isaac/README.md | 86 ++++ .../gno.land/p/wyhaines/rand/isaac/gno.mod | 7 + .../gno.land/p/wyhaines/rand/isaac/isaac.gno | 435 ++++++++++++++++++ .../p/wyhaines/rand/isaac/isaac_test.gno | 165 +++++++ .../p/wyhaines/rand/isaac64/README.md | 97 ++++ .../gno.land/p/wyhaines/rand/isaac64/gno.mod | 7 + .../p/wyhaines/rand/isaac64/isaac64.gno | 429 +++++++++++++++++ .../p/wyhaines/rand/isaac64/isaac64_test.gno | 165 +++++++ .../p/wyhaines/rand/xorshift64star/README.MD | 69 +++ .../p/wyhaines/rand/xorshift64star/gno.mod | 6 + .../rand/xorshift64star/xorshift64star.gno | 172 +++++++ .../xorshift64star/xorshift64star_test.gno | 134 ++++++ .../wyhaines/rand/xorshiftr128plus/README.MD | 60 +++ .../p/wyhaines/rand/xorshiftr128plus/gno.mod | 6 + .../xorshiftr128plus/xorshiftr128plus.gno | 186 ++++++++ .../xorshiftr128plus_test.gno | 142 ++++++ 20 files changed, 2218 insertions(+) create mode 100644 examples/gno.land/p/wyhaines/rand/isaac/README.md create mode 100644 examples/gno.land/p/wyhaines/rand/isaac/gno.mod create mode 100644 examples/gno.land/p/wyhaines/rand/isaac/isaac.gno create mode 100644 examples/gno.land/p/wyhaines/rand/isaac/isaac_test.gno create mode 100644 examples/gno.land/p/wyhaines/rand/isaac64/README.md create mode 100644 examples/gno.land/p/wyhaines/rand/isaac64/gno.mod create mode 100644 examples/gno.land/p/wyhaines/rand/isaac64/isaac64.gno create mode 100644 examples/gno.land/p/wyhaines/rand/isaac64/isaac64_test.gno create mode 100644 examples/gno.land/p/wyhaines/rand/xorshift64star/README.MD create mode 100644 examples/gno.land/p/wyhaines/rand/xorshift64star/gno.mod create mode 100644 examples/gno.land/p/wyhaines/rand/xorshift64star/xorshift64star.gno create mode 100644 examples/gno.land/p/wyhaines/rand/xorshift64star/xorshift64star_test.gno create mode 100644 examples/gno.land/p/wyhaines/rand/xorshiftr128plus/README.MD create mode 100644 examples/gno.land/p/wyhaines/rand/xorshiftr128plus/gno.mod create mode 100644 examples/gno.land/p/wyhaines/rand/xorshiftr128plus/xorshiftr128plus.gno create mode 100644 examples/gno.land/p/wyhaines/rand/xorshiftr128plus/xorshiftr128plus_test.gno diff --git a/examples/gno.land/p/demo/entropy/entropy.gno b/examples/gno.land/p/demo/entropy/entropy.gno index 5e35b8c7227..9e8f656c21b 100644 --- a/examples/gno.land/p/demo/entropy/entropy.gno +++ b/examples/gno.land/p/demo/entropy/entropy.gno @@ -87,3 +87,11 @@ func (i *Instance) Value() uint32 { i.addEntropy() return i.value } + +func (i *Instance) Value64() uint64 { + i.addEntropy() + high := i.value + i.addEntropy() + + return (uint64(high) << 32) | uint64(i.value) +} diff --git a/examples/gno.land/p/demo/entropy/entropy_test.gno b/examples/gno.land/p/demo/entropy/entropy_test.gno index 0deb3ab9aa2..895bfd1e394 100644 --- a/examples/gno.land/p/demo/entropy/entropy_test.gno +++ b/examples/gno.land/p/demo/entropy/entropy_test.gno @@ -33,6 +33,26 @@ func TestInstanceValue(t *testing.T) { } } +func TestInstanceValue64(t *testing.T) { + baseEntropy := New() + baseResult := computeValue64(t, baseEntropy) + + sameHeightEntropy := New() + sameHeightResult := computeValue64(t, sameHeightEntropy) + + if baseResult != sameHeightResult { + t.Errorf("should have the same result: new=%s, base=%s", sameHeightResult, baseResult) + } + + std.TestSkipHeights(1) + differentHeightEntropy := New() + differentHeightResult := computeValue64(t, differentHeightEntropy) + + if baseResult == differentHeightResult { + t.Errorf("should have different result: new=%s, base=%s", differentHeightResult, baseResult) + } +} + func computeValue(t *testing.T, r *Instance) string { t.Helper() @@ -44,3 +64,15 @@ func computeValue(t *testing.T, r *Instance) string { return out } + +func computeValue64(t *testing.T, r *Instance) string { + t.Helper() + + out := "" + for i := 0; i < 10; i++ { + val := int(r.Value64()) + out += strconv.Itoa(val) + " " + } + + return out +} diff --git a/examples/gno.land/p/demo/entropy/z_filetest.gno b/examples/gno.land/p/demo/entropy/z_filetest.gno index 85ed1b10a3d..ddee29b22fd 100644 --- a/examples/gno.land/p/demo/entropy/z_filetest.gno +++ b/examples/gno.land/p/demo/entropy/z_filetest.gno @@ -15,6 +15,7 @@ func main() { println(r.Value()) println(r.Value()) println(r.Value()) + println(r.Value64()) // should be the same println("---") @@ -24,6 +25,7 @@ func main() { println(r.Value()) println(r.Value()) println(r.Value()) + println(r.Value64()) std.TestSkipHeights(1) println("---") @@ -33,6 +35,7 @@ func main() { println(r.Value()) println(r.Value()) println(r.Value()) + println(r.Value64()) } // Output: @@ -42,15 +45,18 @@ func main() { // 1950222777 // 3348280598 // 438354259 +// 6353385488959065197 // --- // 4129293727 // 2141104956 // 1950222777 // 3348280598 // 438354259 +// 6353385488959065197 // --- // 49506731 // 1539580078 // 2695928529 // 1895482388 // 3462727799 +// 16745038698684748445 diff --git a/examples/gno.land/p/demo/ufmt/ufmt_test.gno b/examples/gno.land/p/demo/ufmt/ufmt_test.gno index 1cb7231a611..1a4d4e7e6f2 100644 --- a/examples/gno.land/p/demo/ufmt/ufmt_test.gno +++ b/examples/gno.land/p/demo/ufmt/ufmt_test.gno @@ -44,6 +44,12 @@ func TestSprintf(t *testing.T) { {"uint32 [%v]", []interface{}{uint32(32)}, "uint32 [32]"}, {"uint64 [%d]", []interface{}{uint64(64)}, "uint64 [64]"}, {"uint64 [%v]", []interface{}{uint64(64)}, "uint64 [64]"}, + {"float64 [%e]", []interface{}{float64(64.1)}, "float64 [6.41e+01]"}, + {"float64 [%E]", []interface{}{float64(64.1)}, "float64 [6.41E+01]"}, + {"float64 [%f]", []interface{}{float64(64.1)}, "float64 [64.100000]"}, + {"float64 [%F]", []interface{}{float64(64.1)}, "float64 [64.100000]"}, + {"float64 [%g]", []interface{}{float64(64.1)}, "float64 [64.1]"}, + {"float64 [%G]", []interface{}{float64(64.1)}, "float64 [64.1]"}, {"bool [%t]", []interface{}{true}, "bool [true]"}, {"bool [%v]", []interface{}{true}, "bool [true]"}, {"bool [%t]", []interface{}{false}, "bool [false]"}, diff --git a/examples/gno.land/p/wyhaines/rand/isaac/README.md b/examples/gno.land/p/wyhaines/rand/isaac/README.md new file mode 100644 index 00000000000..05f4a94425f --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/isaac/README.md @@ -0,0 +1,86 @@ +# package isaac // import "gno.land/p/demo/math/rand/isaac" + +This is a port of the ISAAC cryptographically secure PRNG, +originally based on the reference implementation found at +https://burtleburtle.net/bob/rand/isaacafa.html + +ISAAC has excellent statistical properties, with long cycle times, and +uniformly distributed, unbiased, and unpredictable number generation. It can +not be distinguished from real random data, and in three decades of scrutiny, +no practical attacks have been found. + +The default random number algorithm in gno was ported from Go's v2 rand +implementatoon, which defaults to the PCG algorithm. This algorithm is +commonly used in language PRNG implementations because it has modest seeding +requirements, and generates statistically strong randomness. + +This package provides an implementation of the 32-bit ISAAC PRNG algorithm. This +algorithm provides very strong statistical performance, and is cryptographically +secure, while still being substantially faster than the default PCG +implementation in `math/rand`. Note that this package does implement a `Uint64()` +function in order to generate a 64 bit number out of two 32 bit numbers. Doing this +makes the generator only slightly faster than PCG, however, + +Note that the approach to seeing with ISAAC is very important for best results, +and seeding with ISAAC is not as simple as seeding with a single uint64 value. +The ISAAC algorithm requires a 256-element seed. If used for cryptographic +purposes, this will likely require entropy generated off-chain for actual +cryptographically secure seeding. For other purposes, however, one can utilize +the built-in seeding mechanism, which will leverage the xorshiftr128plus PRNG to +generate any missing seeds if fewer than 256 are provided. + + +``` +Benchmark +--------- +PCG: 1000000 Uint64 generated in 15.58s +ISAAC: 1000000 Uint64 generated in 13.23s (uint64) +ISAAC: 1000000 Uint32 generated in 6.43s (uint32) +Ratio: x1.18 times faster than PCG (uint64) +Ratio: x2.42 times faster than PCG (uint32) +``` + +Use it directly: + +``` +prng = isaac.New() // pass 0 to 256 uint32 seeds; if fewer than 256 are provided, the rest + // will be generated using the xorshiftr128plus PRNG. +``` + +Or use it as a drop-in replacement for the default PRNT in Rand: + +``` +source = isaac.New() +prng := rand.New(source) +``` + +# TYPES + +` +type ISAAC struct { + // Has unexported fields. +} +` + +`func New(seeds ...uint32) *ISAAC` + ISAAC requires a large, 256-element seed. This implementation will leverage + the entropy package combined with the the xorshiftr128plus PRNG to generate + any missing seeds of fewer than the required number of arguments are + provided. + +`func (isaac *ISAAC) MarshalBinary() ([]byte, error)` + MarshalBinary() returns a byte array that encodes the state of the PRNG. + This can later be used with UnmarshalBinary() to restore the state of the + PRNG. MarshalBinary implements the encoding.BinaryMarshaler interface. + +`func (isaac *ISAAC) Seed(seed [256]uint32)` + +`func (isaac *ISAAC) Uint32() uint32` + +`func (isaac *ISAAC) Uint64() uint64` + +`func (isaac *ISAAC) UnmarshalBinary(data []byte) error` + UnmarshalBinary() restores the state of the PRNG from a byte array + that was created with MarshalBinary(). UnmarshalBinary implements the + encoding.BinaryUnmarshaler interface. + diff --git a/examples/gno.land/p/wyhaines/rand/isaac/gno.mod b/examples/gno.land/p/wyhaines/rand/isaac/gno.mod new file mode 100644 index 00000000000..0cca6aa5174 --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/isaac/gno.mod @@ -0,0 +1,7 @@ +module gno.land/p/wyhaines/rand/isaac + +require ( + gno.land/p/demo/entropy v0.0.0-latest + gno.land/p/demo/ufmt v0.0.0-latest + gno.land/p/wyhaines/rand/xorshiftr128plus v0.0.0-latest +) diff --git a/examples/gno.land/p/wyhaines/rand/isaac/isaac.gno b/examples/gno.land/p/wyhaines/rand/isaac/isaac.gno new file mode 100644 index 00000000000..4508dd5d5af --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/isaac/isaac.gno @@ -0,0 +1,435 @@ +// This is a port of the ISAAC cryptographically secure PRNG, originally based on the reference +// implementation found at https://burtleburtle.net/bob/rand/isaacafa.html +// +// ISAAC has excellent statistical properties, with long cycle times, and uniformly distributed, +// unbiased, and unpredictable number generation. It can not be distinguished from real random +// data, and in three decades of scrutiny, no practical attacks have been found. +// +// The default random number algorithm in gno was ported from Go's v2 rand implementation, which +// defaults to the PCG algorithm. This algorithm is commonly used in language PRNG implementations +// because it has modest seeding requirements, and generates statistically strong randomness. +// +// This package provides an implementation of the 32-bit ISAAC PRNG algorithm. This +// algorithm provides very strong statistical performance, and is cryptographically +// secure, while still being substantially faster than the default PCG +// implementation in `math/rand`. Note that this package does implement a `Uint64()` +// function in order to generate a 64 bit number out of two 32 bit numbers. Doing this +// makes the generator only slightly faster than PCG, however, +// +// Note that the approach to seeing with ISAAC is very important for best results, and seeding with +// ISAAC is not as simple as seeding with a single uint64 value. The ISAAC algorithm requires a +// 256-element seed. If used for cryptographic purposes, this will likely require entropy generated +// off-chain for actual cryptographically secure seeding. For other purposes, however, one can +// utilize the built-in seeding mechanism, which will leverage the xorshiftr128plus PRNG to generate +// any missing seeds if fewer than 256 are provided. +// +// Benchmark +// --------- +// PCG: 1000000 Uint64 generated in 15.58s +// ISAAC: 1000000 Uint64 generated in 13.23s +// ISAAC: 1000000 Uint32 generated in 6.43s +// Ratio: x1.18 times faster than PCG (uint64) +// Ratio: x2.42 times faster than PCG (uint32) +// +// Use it directly: +// +// prng = isaac.New() // pass 0 to 256 uint32 seeds; if fewer than 256 are provided, the rest +// // will be generated using the xorshiftr128plus PRNG. +// +// Or use it as a drop-in replacement for the default PRNG in Rand: +// +// source = isaac.New() +// prng := rand.New(source) +package isaac + +import ( + "errors" + "math" + "math/rand" + + "gno.land/p/demo/entropy" + "gno.land/p/demo/ufmt" + "gno.land/p/wyhaines/rand/xorshiftr128plus" +) + +type ISAAC struct { + randrsl [256]uint32 + randcnt uint32 + mm [256]uint32 + aa, bb, cc uint32 + seed [256]uint32 +} + +// ISAAC requires a large, 256-element seed. This implementation will leverage the entropy +// package combined with the the xorshiftr128plus PRNG to generate any missing seeds of +// fewer than the required number of arguments are provided. +func New(seeds ...uint32) *ISAAC { + isaac := &ISAAC{} + seed := [256]uint32{} + + index := 0 + for index = 0; index < len(seeds); index++ { + seed[index] = seeds[index] + } + + if index < 4 { + e := entropy.New() + for ; index < 4; index++ { + seed[index] = e.Value() + } + } + + // Use up to the first four seeds as seeding inputs for xorshiftr128+, in order to + // use it to provide any remaining missing seeds. + prng := xorshiftr128plus.New( + (uint64(seed[0])<<32)|uint64(seed[1]), + (uint64(seed[2])<<32)|uint64(seed[3]), + ) + for ; index < 256; index += 2 { + val := prng.Uint64() + seed[index] = uint32(val & 0xffffffff) + if index+1 < 256 { + seed[index+1] = uint32(val >> 32) + } + } + isaac.Seed(seed) + return isaac +} + +func (isaac *ISAAC) Seed(seed [256]uint32) { + isaac.randrsl = seed + isaac.seed = seed + isaac.randinit(true) +} + +// beUint32() decodes a uint32 from a set of four bytes, assuming big endian encoding. +// binary.bigEndian.Uint32, copied to avoid dependency +func beUint32(b []byte) uint32 { + _ = b[3] // bounds check hint to compiler; see golang.org/issue/14808 + return uint32(b[3]) | uint32(b[2])<<8 | uint32(b[1])<<16 | uint32(b[0])<<24 +} + +// bePutUint32() encodes a uint64 into a buffer of eight bytes. +// binary.bigEndian.PutUint32, copied to avoid dependency +func bePutUint32(b []byte, v uint32) { + _ = b[3] // early bounds check to guarantee safety of writes below + b[0] = byte(v >> 24) + b[1] = byte(v >> 16) + b[2] = byte(v >> 8) + b[3] = byte(v) +} + +// A label to identify the marshalled data. +var marshalISAACLabel = []byte("isaac:") + +// MarshalBinary() returns a byte array that encodes the state of the PRNG. This can later be used +// with UnmarshalBinary() to restore the state of the PRNG. +// MarshalBinary implements the encoding.BinaryMarshaler interface. +func (isaac *ISAAC) MarshalBinary() ([]byte, error) { + b := make([]byte, 3094) // 6 + 1024 + 1024 + 1024 + 4 + 4 + 4 + 4 == 3090 + copy(b, marshalISAACLabel) + for i := 0; i < 256; i++ { + bePutUint32(b[6+i*4:], isaac.seed[i]) + } + for i := 256; i < 512; i++ { + bePutUint32(b[6+i*4:], isaac.randrsl[i-256]) + } + for i := 512; i < 768; i++ { + bePutUint32(b[6+i*4:], isaac.mm[i-512]) + } + bePutUint32(b[3078:], isaac.aa) + bePutUint32(b[3082:], isaac.bb) + bePutUint32(b[3086:], isaac.cc) + bePutUint32(b[3090:], isaac.randcnt) + + return b, nil +} + +// errUnmarshalISAAC is returned when unmarshalling fails. +var errUnmarshalISAAC = errors.New("invalid ISAAC encoding") + +// UnmarshalBinary() restores the state of the PRNG from a byte array that was created with MarshalBinary(). +// UnmarshalBinary implements the encoding.BinaryUnmarshaler interface. +func (isaac *ISAAC) UnmarshalBinary(data []byte) error { + if len(data) != 3094 || string(data[:6]) != string(marshalISAACLabel) { + return errUnmarshalISAAC + } + for i := 0; i < 256; i++ { + isaac.seed[i] = beUint32(data[6+i*4:]) + } + for i := 256; i < 512; i++ { + isaac.randrsl[i-256] = beUint32(data[6+i*4:]) + } + for i := 512; i < 768; i++ { + isaac.mm[i-512] = beUint32(data[6+i*4:]) + } + isaac.aa = beUint32(data[3078:]) + isaac.bb = beUint32(data[3082:]) + isaac.cc = beUint32(data[3086:]) + isaac.randcnt = beUint32(data[3090:]) + return nil +} + +func (isaac *ISAAC) randinit(flag bool) { + isaac.aa = 0 + isaac.bb = 0 + isaac.cc = 0 + + var a, b, c, d, e, f, g, h uint32 = 0x9e3779b9, 0x9e3779b9, 0x9e3779b9, 0x9e3779b9, 0x9e3779b9, 0x9e3779b9, 0x9e3779b9, 0x9e3779b9 + + for i := 0; i < 4; i++ { + a ^= b << 11 + d += a + b += c + b ^= c >> 2 + e += b + c += d + c ^= d << 8 + f += c + d += e + d ^= e >> 16 + g += d + e += f + e ^= f << 10 + h += e + f += g + f ^= g >> 4 + a += f + g += h + g ^= h << 8 + b += g + h += a + h ^= a >> 9 + c += h + a += b + } + + for i := 0; i < 256; i += 8 { + if flag { + a += isaac.randrsl[i] + b += isaac.randrsl[i+1] + c += isaac.randrsl[i+2] + d += isaac.randrsl[i+3] + e += isaac.randrsl[i+4] + f += isaac.randrsl[i+5] + g += isaac.randrsl[i+6] + h += isaac.randrsl[i+7] + } + + a ^= b << 11 + d += a + b += c + b ^= c >> 2 + e += b + c += d + c ^= d << 8 + f += c + d += e + d ^= e >> 16 + g += d + e += f + e ^= f << 10 + h += e + f += g + f ^= g >> 4 + a += f + g += h + g ^= h << 8 + b += g + h += a + h ^= a >> 9 + c += h + a += b + + isaac.mm[i] = a + isaac.mm[i+1] = b + isaac.mm[i+2] = c + isaac.mm[i+3] = d + isaac.mm[i+4] = e + isaac.mm[i+5] = f + isaac.mm[i+6] = g + isaac.mm[i+7] = h + } + + if flag { + for i := 0; i < 256; i += 8 { + a += isaac.mm[i] + b += isaac.mm[i+1] + c += isaac.mm[i+2] + d += isaac.mm[i+3] + e += isaac.mm[i+4] + f += isaac.mm[i+5] + g += isaac.mm[i+6] + h += isaac.mm[i+7] + + a ^= b << 11 + d += a + b += c + b ^= c >> 2 + e += b + c += d + c ^= d << 8 + f += c + d += e + d ^= e >> 16 + g += d + e += f + e ^= f << 10 + h += e + f += g + f ^= g >> 4 + a += f + g += h + g ^= h << 8 + b += g + h += a + h ^= a >> 9 + c += h + a += b + + isaac.mm[i] = a + isaac.mm[i+1] = b + isaac.mm[i+2] = c + isaac.mm[i+3] = d + isaac.mm[i+4] = e + isaac.mm[i+5] = f + isaac.mm[i+6] = g + isaac.mm[i+7] = h + } + } + + isaac.isaac() + isaac.randcnt = uint32(256) +} + +func (isaac *ISAAC) isaac() { + isaac.cc++ + isaac.bb += isaac.cc + + for i := 0; i < 256; i++ { + x := isaac.mm[i] + switch i % 4 { + case 0: + isaac.aa ^= isaac.aa << 13 + case 1: + isaac.aa ^= isaac.aa >> 6 + case 2: + isaac.aa ^= isaac.aa << 2 + case 3: + isaac.aa ^= isaac.aa >> 16 + } + isaac.aa += isaac.mm[(i+128)&0xff] + + y := isaac.mm[(x>>2)&0xff] + isaac.aa + isaac.bb + isaac.mm[i] = y + isaac.bb = isaac.mm[(y>>10)&0xff] + x + isaac.randrsl[i] = isaac.bb + } +} + +// Returns a random uint32. +func (isaac *ISAAC) Uint32() uint32 { + if isaac.randcnt == uint32(0) { + isaac.isaac() + isaac.randcnt = uint32(256) + } + isaac.randcnt-- + return isaac.randrsl[isaac.randcnt] +} + +// Returns a random uint64 by combining two uint32s. +func (isaac *ISAAC) Uint64() uint64 { + return uint64(isaac.Uint32()) | (uint64(isaac.Uint32()) << 32) +} + +// Until there is better benchmarking support in gno, you can test the performance of this PRNG with this function. +// This isn't perfect, since it will include the startup time of gno in the results, but this will give you a timing +// for generating a million random uint64 numbers on any unix based system: +// +// `time gno run -expr 'benchmarkISAAC()' xorshift64star.gno +func benchmarkISAAC(_iterations ...int) { + iterations := 1000000 + if len(_iterations) > 0 { + iterations = _iterations[0] + } + isaac := New() + + for i := 0; i < iterations; i++ { + _ = isaac.Uint64() + } + ufmt.Println(ufmt.Sprintf("ISAAC: generate %d uint64\n", iterations)) +} + +// The averageISAAC() function is a simple benchmarking helper to demonstrate +// the most basic statistical property of the ISAAC PRNG. +func averageISAAC(_iterations ...int) { + target := uint64(500000) + iterations := 1000000 + var squares [1000000]uint64 + + ufmt.Println( + ufmt.Sprintf( + "Averaging %d random numbers. The average should be very close to %d.\n", + iterations, + target)) + + if len(_iterations) > 0 { + iterations = _iterations[0] + } + isaac := New(987654321, 123456789, 999999999, 111111111) + + var average float64 = 0 + for i := 0; i < iterations; i++ { + n := isaac.Uint64()%(target*2) + 1 + average += (float64(n) - average) / float64(i+1) + squares[i] = n + } + + sum_of_squares := uint64(0) + // transform numbers into their squares of the distance from the average + for i := 0; i < iterations; i++ { + difference := average - float64(squares[i]) + square := uint64(difference * difference) + sum_of_squares += square + } + + ufmt.Println(ufmt.Sprintf("ISAAC average of %d uint64: %f\n", iterations, average)) + ufmt.Println(ufmt.Sprintf("ISAAC standard deviation : %f\n", math.Sqrt(float64(sum_of_squares)/float64(iterations)))) + ufmt.Println(ufmt.Sprintf("ISAAC theoretical perfect deviation: %f\n", (float64(target*2)-1)/math.Sqrt(12))) +} + +func averagePCG(_iterations ...int) { + target := uint64(500000) + iterations := 1000000 + var squares [1000000]uint64 + + ufmt.Println( + ufmt.Sprintf( + "Averaging %d random numbers. The average should be very close to %d.\n", + iterations, + target)) + + if len(_iterations) > 0 { + iterations = _iterations[0] + } + isaac := rand.NewPCG(987654321, 123456789) + + var average float64 = 0 + for i := 0; i < iterations; i++ { + n := isaac.Uint64()%(target*2) + 1 + average += (float64(n) - average) / float64(i+1) + squares[i] = n + } + + sum_of_squares := uint64(0) + // transform numbers into their squares of the distance from the average + for i := 0; i < iterations; i++ { + difference := average - float64(squares[i]) + square := uint64(difference * difference) + sum_of_squares += square + } + + ufmt.Println(ufmt.Sprintf("PCG average of %d uint64: %f\n", iterations, average)) + ufmt.Println(ufmt.Sprintf("PCG standard deviation : %f\n", math.Sqrt(float64(sum_of_squares)/float64(iterations)))) + ufmt.Println(ufmt.Sprintf("PCG theoretical perfect deviation: %f\n", (float64(target*2)-1)/math.Sqrt(12))) +} diff --git a/examples/gno.land/p/wyhaines/rand/isaac/isaac_test.gno b/examples/gno.land/p/wyhaines/rand/isaac/isaac_test.gno new file mode 100644 index 00000000000..b08621e271c --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/isaac/isaac_test.gno @@ -0,0 +1,165 @@ +package isaac + +import ( + "math/rand" + "testing" +) + +type OpenISAAC struct { + Randrsl [256]uint32 + Randcnt uint32 + Mm [256]uint32 + Aa, Bb, Cc uint32 + Seed [256]uint32 +} + +func TestISAACSeeding(t *testing.T) { + isaac := New() +} + +func TestISAACRand(t *testing.T) { + source := New(987654321) + rng := rand.New(source) + + // Expected outputs for the first 5 random floats with the given seed + expected := []float64{ + 0.17828173023837635, + 0.7327795780287832, + 0.4850369074875177, + 0.9474842397428482, + 0.6747135561813891, + 0.7522507082868403, + 0.041115261836534356, + 0.7405243709084567, + 0.672863376128768, + 0.11866211399980553, + } + + for i, exp := range expected { + val := rng.Float64() + if exp != val { + t.Errorf("Rand.Float64() at iteration %d: got %g, expected %g", i, val, exp) + } + } +} + +func TestISAACUint64(t *testing.T) { + isaac := New() + + expected := []uint64{ + 5986068031949215749, + 10437354066128700566, + 13478007513323023970, + 8969511410255984224, + 3869229557962857982, + 1762449743873204415, + 5292356290662282456, + 7893982194485405616, + 4296136494566588699, + 12414349056998262772, + } + + for i, exp := range expected { + val := isaac.Uint64() + if exp != val { + t.Errorf("ISAAC.Uint64() at iteration %d: got %d, expected %d", i, val, exp) + } + } +} + +func dupState(i *ISAAC) *OpenISAAC { + state := &OpenISAAC{} + state.Seed = i.seed + state.Randrsl = i.randrsl + state.Mm = i.mm + state.Aa = i.aa + state.Bb = i.bb + state.Cc = i.cc + state.Randcnt = i.randcnt + + return state +} + +func TestISAACMarshalUnmarshal(t *testing.T) { + isaac := New() + + expected1 := []uint64{ + 5986068031949215749, + 10437354066128700566, + 13478007513323023970, + 8969511410255984224, + 3869229557962857982, + } + + expected2 := []uint64{ + 1762449743873204415, + 5292356290662282456, + 7893982194485405616, + 4296136494566588699, + 12414349056998262772, + } + + for i, exp := range expected1 { + val := isaac.Uint64() + if exp != val { + t.Errorf("ISAAC.Uint64() at iteration %d: got %d, expected %d", i, val, exp) + } + } + + marshalled, err := isaac.MarshalBinary() + + t.Logf("State: [%v]\n", dupState(isaac)) + t.Logf("Marshalled State: [%x] -- %v\n", marshalled, err) + state_before := dupState(isaac) + + if err != nil { + t.Errorf("ISAAC.MarshalBinary() error: %v", err) + } + + // Advance state by one number; then check the next 5. The expectation is that they _will_ fail. + isaac.Uint64() + + for i, exp := range expected2 { + val := isaac.Uint64() + if exp == val { + t.Errorf(" Iteration %d matched %d; which is from iteration %d; something strange is happening.", (i + 6), val, (i + 5)) + } + } + + t.Logf("State before unmarshall: [%v]\n", dupState(isaac)) + + // Now restore the state of the PRNG + err = isaac.UnmarshalBinary(marshalled) + + t.Logf("State after unmarshall: [%v]\n", dupState(isaac)) + + if state_before.Seed != dupState(isaac).Seed { + t.Errorf("Seed mismatch") + } + if state_before.Randrsl != dupState(isaac).Randrsl { + t.Errorf("Randrsl mismatch") + } + if state_before.Mm != dupState(isaac).Mm { + t.Errorf("Mm mismatch") + } + if state_before.Aa != dupState(isaac).Aa { + t.Errorf("Aa mismatch") + } + if state_before.Bb != dupState(isaac).Bb { + t.Errorf("Bb mismatch") + } + if state_before.Cc != dupState(isaac).Cc { + t.Errorf("Cc mismatch") + } + if state_before.Randcnt != dupState(isaac).Randcnt { + t.Errorf("Randcnt mismatch") + } + + // Now we should be back on track for the last 5 numbers + for i, exp := range expected2 { + val := isaac.Uint64() + if exp != val { + t.Errorf("ISAAC.Uint64() at iteration %d: got %d, expected %d", (i + 5), val, exp) + } + } +} diff --git a/examples/gno.land/p/wyhaines/rand/isaac64/README.md b/examples/gno.land/p/wyhaines/rand/isaac64/README.md new file mode 100644 index 00000000000..813b062a5cd --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/isaac64/README.md @@ -0,0 +1,97 @@ +# package isaac64 // import "gno.land/p/demo/math/rand/isaac64" + +This is a port of the 64-bit version of the ISAAC cryptographically +secure PRNG, originally based on the reference implementation found at +https://burtleburtle.net/bob/rand/isaacafa.html + +ISAAC has excellent statistical properties, with long cycle times, and +uniformly distributed, unbiased, and unpredictable number generation. It can +not be distinguished from real random data, and in three decades of scrutiny, +no practical attacks have been found. + +The default random number algorithm in gno was ported from Go's v2 rand +implementatoon, which defaults to the PCG algorithm. This algorithm is +commonly used in language PRNG implementations because it has modest seeding +requirements, and generates statistically strong randomness. + +This package provides an implementation of the 64-bit ISAAC PRNG algorithm. This +algorithm provides very strong statistical performance, and is cryptographically +secure, while still being substantially faster than the default PCG +implementation in `math/rand`. + +Note that the approach to seeing with ISAAC is very important for best results, +and seeding with ISAAC is not as simple as seeding with a single uint64 value. +The ISAAC algorithm requires a 256-element seed. If used for cryptographic +purposes, this will likely require entropy generated off-chain for actual +cryptographically secure seeding. For other purposes, however, one can utilize +the built-in seeding mechanism, which will leverage the xorshiftr128plus PRNG to +generate any missing seeds if fewer than 256 are provided. + + +``` +Benchmark +--------- +PCG: 1000000 Uint64 generated in 15.58s +ISAAC: 1000000 Uint64 generated in 8.95s +ISAAC: 1000000 Uint32 generated in 7.66s +Ratio: x1.74 times faster than PCG (uint64) +Ratio: x2.03 times faster than PCG (uint32) +``` + +Use it directly: + + +``` +prng = isaac.New() // pass 0 to 256 uint64 seeds; if fewer than 256 are provided, the rest + // will be generated using the xorshiftr128plus PRNG. +``` + +Or use it as a drop-in replacement for the default PRNT in Rand: + +``` +source = isaac64.New() +prng := rand.New(source) +``` + +## CONSTANTS + + +``` +const ( + RANDSIZL = 8 + RANDSIZ = 1 << RANDSIZL // 256 +) +``` + +## TYPES + + +``` +type ISAAC struct { + // Has unexported fields. +} +``` + +`func New(seeds ...uint64) *ISAAC` +ISAAC requires a large, 256-element seed. This implementation will leverage +the entropy package combined with the xorshiftr128plus PRNG to generate any +missing seeds if fewer than the required number of arguments are provided. + +`func (isaac *ISAAC) MarshalBinary() ([]byte, error)` +MarshalBinary() returns a byte array that encodes the state of the PRNG. +This can later be used with UnmarshalBinary() to restore the state of the +PRNG. MarshalBinary implements the encoding.BinaryMarshaler interface. + +`func (isaac *ISAAC) Seed(seed [256]uint64)` +Reinitialize the generator with a new seed. A seed must be composed of 256 uint64. + +`func (isaac *ISAAC) Uint32() uint32` +Return a 32 bit random integer, composed of the high 32 bits of the generated 32 bit result. + +`func (isaac *ISAAC) Uint64() uint64` +Return a 64 bit random integer. + +`func (isaac *ISAAC) UnmarshalBinary(data []byte) error` +UnmarshalBinary() restores the state of the PRNG from a byte array +that was created with MarshalBinary(). UnmarshalBinary implements the +encoding.BinaryUnmarshaler interface. diff --git a/examples/gno.land/p/wyhaines/rand/isaac64/gno.mod b/examples/gno.land/p/wyhaines/rand/isaac64/gno.mod new file mode 100644 index 00000000000..dbc8713094e --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/isaac64/gno.mod @@ -0,0 +1,7 @@ +module gno.land/p/wyhaines/rand/isaac64 + +require ( + gno.land/p/demo/entropy v0.0.0-latest + gno.land/p/demo/ufmt v0.0.0-latest + gno.land/p/wyhaines/rand/xorshiftr128plus v0.0.0-latest +) diff --git a/examples/gno.land/p/wyhaines/rand/isaac64/isaac64.gno b/examples/gno.land/p/wyhaines/rand/isaac64/isaac64.gno new file mode 100644 index 00000000000..6f2d95150fc --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/isaac64/isaac64.gno @@ -0,0 +1,429 @@ +// This is a port of the 64-bit version of the ISAAC cryptographically secure PRNG, originally +// based on the reference implementation found at https://burtleburtle.net/bob/rand/isaacafa.html +// +// ISAAC has excellent statistical properties, with long cycle times, and uniformly distributed, +// unbiased, and unpredictable number generation. It can not be distinguished from real random +// data, and in three decades of scrutiny, no practical attacks have been found. +// +// The default random number algorithm in gno was ported from Go's v2 rand implementatoon, which +// defaults to the PCG algorithm. This algorithm is commonly used in language PRNG implementations +// because it has modest seeding requirements, and generates statistically strong randomness. +// +// This package provides an implementation of the 64-bit ISAAC PRNG algorithm. This algorithm +// provides very strong statistical performance, and is cryptographically secure, while still +// being substantially faster than the default PCG implementation in `math/rand`. +// +// Note that the approach to seeing with ISAAC is very important for best results, and seeding with +// ISAAC is not as simple as seeding with a single uint64 value. The ISAAC algorithm requires a +// 256-element seed. If used for cryptographic purposes, this will likely require entropy generated +// off-chain for actual cryptographically secure seeding. For other purposes, however, one can +// utilize the built-in seeding mechanism, which will leverage the xorshiftr128plus PRNG to generate +// any missing seeds if fewer than 256 are provided. +// +// Benchmark +// --------- +// PCG: 1000000 Uint64 generated in 15.58s +// ISAAC: 1000000 Uint64 generated in 8.95s +// ISAAC: 1000000 Uint32 generated in 7.66s +// Ratio: x1.74 times faster than PCG (uint64) +// Ratio: x2.03 times faster than PCG (uint32) +// +// Use it directly: +// +// prng = isaac.New() // pass 0 to 256 uint64 seeds; if fewer than 256 are provided, the rest +// // will be generated using the xorshiftr128plus PRNG. +// +// Or use it as a drop-in replacement for the default PRNT in Rand: +// +// source = isaac64.New() +// prng := rand.New(source) +package isaac64 + +import ( + "errors" + "math" + + "gno.land/p/demo/entropy" + "gno.land/p/demo/ufmt" + "gno.land/p/wyhaines/rand/xorshiftr128plus" +) + +const ( + RANDSIZL = 8 + RANDSIZ = 1 << RANDSIZL // 256 +) + +type ISAAC struct { + randrsl [256]uint64 + randcnt uint64 + mm [256]uint64 + aa, bb, cc uint64 + seed [256]uint64 +} + +// ISAAC requires a large, 256-element seed. This implementation will leverage the entropy +// package combined with the xorshiftr128plus PRNG to generate any missing seeds if fewer than +// the required number of arguments are provided. +func New(seeds ...uint64) *ISAAC { + isaac := &ISAAC{} + seed := [256]uint64{} + + index := 0 + for index = 0; index < len(seeds) && index < 256; index++ { + seed[index] = seeds[index] + } + + if index < 2 { + e := entropy.New() + for ; index < 2; index++ { + seed[index] = e.Value64() + } + } + + // Use the first two seeds as seeding inputs for xorshiftr128plus, in order to + // use it to provide any remaining missing seeds. + prng := xorshiftr128plus.New( + seed[0], + seed[1], + ) + for ; index < 256; index++ { + seed[index] = prng.Uint64() + } + isaac.Seed(seed) + return isaac +} + +// Reinitialize the generator with a new seed. A seed must be composed of 256 uint64. +func (isaac *ISAAC) Seed(seed [256]uint64) { + isaac.randrsl = seed + isaac.seed = seed + isaac.randinit(true) +} + +// beUint64() decodes a uint64 from a set of eight bytes, assuming big endian encoding. +func beUint64(b []byte) uint64 { + _ = b[7] // bounds check hint to compiler + return uint64(b[7]) | uint64(b[6])<<8 | uint64(b[5])<<16 | uint64(b[4])<<24 | + uint64(b[3])<<32 | uint64(b[2])<<40 | uint64(b[1])<<48 | uint64(b[0])<<56 +} + +// bePutUint64() encodes a uint64 into a buffer of eight bytes. +func bePutUint64(b []byte, v uint64) { + _ = b[7] // early bounds check to guarantee safety of writes below + b[0] = byte(v >> 56) + b[1] = byte(v >> 48) + b[2] = byte(v >> 40) + b[3] = byte(v >> 32) + b[4] = byte(v >> 24) + b[5] = byte(v >> 16) + b[6] = byte(v >> 8) + b[7] = byte(v) +} + +// A label to identify the marshalled data. +var marshalISAACLabel = []byte("isaac:") + +// MarshalBinary() returns a byte array that encodes the state of the PRNG. This can later be used +// with UnmarshalBinary() to restore the state of the PRNG. +// MarshalBinary implements the encoding.BinaryMarshaler interface. +func (isaac *ISAAC) MarshalBinary() ([]byte, error) { + b := make([]byte, 6+2048*3+8*3+8) // 6 + 2048*3 + 8*3 + 8 == 6182 + copy(b, marshalISAACLabel) + offset := 6 + for i := 0; i < 256; i++ { + bePutUint64(b[offset:], isaac.seed[i]) + offset += 8 + } + for i := 0; i < 256; i++ { + bePutUint64(b[offset:], isaac.randrsl[i]) + offset += 8 + } + for i := 0; i < 256; i++ { + bePutUint64(b[offset:], isaac.mm[i]) + offset += 8 + } + bePutUint64(b[offset:], isaac.aa) + offset += 8 + bePutUint64(b[offset:], isaac.bb) + offset += 8 + bePutUint64(b[offset:], isaac.cc) + offset += 8 + bePutUint64(b[offset:], isaac.randcnt) + return b, nil +} + +// errUnmarshalISAAC is returned when unmarshalling fails. +var errUnmarshalISAAC = errors.New("invalid ISAAC encoding") + +// UnmarshalBinary() restores the state of the PRNG from a byte array that was created with MarshalBinary(). +// UnmarshalBinary implements the encoding.BinaryUnmarshaler interface. +func (isaac *ISAAC) UnmarshalBinary(data []byte) error { + if len(data) != 6182 || string(data[:6]) != string(marshalISAACLabel) { + return errUnmarshalISAAC + } + offset := 6 + for i := 0; i < 256; i++ { + isaac.seed[i] = beUint64(data[offset:]) + offset += 8 + } + for i := 0; i < 256; i++ { + isaac.randrsl[i] = beUint64(data[offset:]) + offset += 8 + } + for i := 0; i < 256; i++ { + isaac.mm[i] = beUint64(data[offset:]) + offset += 8 + } + isaac.aa = beUint64(data[offset:]) + offset += 8 + isaac.bb = beUint64(data[offset:]) + offset += 8 + isaac.cc = beUint64(data[offset:]) + offset += 8 + isaac.randcnt = beUint64(data[offset:]) + return nil +} + +func (isaac *ISAAC) randinit(flag bool) { + var a, b, c, d, e, f, g, h uint64 + isaac.aa = 0 + isaac.bb = 0 + isaac.cc = 0 + + a = 0x9e3779b97f4a7c13 + b = 0x9e3779b97f4a7c13 + c = 0x9e3779b97f4a7c13 + d = 0x9e3779b97f4a7c13 + e = 0x9e3779b97f4a7c13 + f = 0x9e3779b97f4a7c13 + g = 0x9e3779b97f4a7c13 + h = 0x9e3779b97f4a7c13 + + // scramble it + for i := 0; i < 4; i++ { + mix(&a, &b, &c, &d, &e, &f, &g, &h) + } + + // fill in mm[] with messy stuff + for i := 0; i < RANDSIZ; i += 8 { + if flag { + a += isaac.randrsl[i] + b += isaac.randrsl[i+1] + c += isaac.randrsl[i+2] + d += isaac.randrsl[i+3] + e += isaac.randrsl[i+4] + f += isaac.randrsl[i+5] + g += isaac.randrsl[i+6] + h += isaac.randrsl[i+7] + } + mix(&a, &b, &c, &d, &e, &f, &g, &h) + isaac.mm[i] = a + isaac.mm[i+1] = b + isaac.mm[i+2] = c + isaac.mm[i+3] = d + isaac.mm[i+4] = e + isaac.mm[i+5] = f + isaac.mm[i+6] = g + isaac.mm[i+7] = h + } + + if flag { + // do a second pass to make all of the seed affect all of mm + for i := 0; i < RANDSIZ; i += 8 { + a += isaac.mm[i] + b += isaac.mm[i+1] + c += isaac.mm[i+2] + d += isaac.mm[i+3] + e += isaac.mm[i+4] + f += isaac.mm[i+5] + g += isaac.mm[i+6] + h += isaac.mm[i+7] + mix(&a, &b, &c, &d, &e, &f, &g, &h) + isaac.mm[i] = a + isaac.mm[i+1] = b + isaac.mm[i+2] = c + isaac.mm[i+3] = d + isaac.mm[i+4] = e + isaac.mm[i+5] = f + isaac.mm[i+6] = g + isaac.mm[i+7] = h + } + } + + isaac.isaac() + isaac.randcnt = RANDSIZ +} + +func mix(a, b, c, d, e, f, g, h *uint64) { + *a -= *e + *f ^= *h >> 9 + *h += *a + + *b -= *f + *g ^= *a << 9 + *a += *b + + *c -= *g + *h ^= *b >> 23 + *b += *c + + *d -= *h + *a ^= *c << 15 + *c += *d + + *e -= *a + *b ^= *d >> 14 + *d += *e + + *f -= *b + *c ^= *e << 20 + *e += *f + + *g -= *c + *d ^= *f >> 17 + *f += *g + + *h -= *d + *e ^= *g << 14 + *g += *h +} + +func ind(mm []uint64, x uint64) uint64 { + return mm[(x>>3)&(RANDSIZ-1)] +} + +func (isaac *ISAAC) isaac() { + var a, b, x, y uint64 + a = isaac.aa + b = isaac.bb + isaac.cc + 1 + isaac.cc++ + + m := isaac.mm[:] + r := isaac.randrsl[:] + + var i, m2Index int + + // First half + for i = 0; i < RANDSIZ/2; i++ { + m2Index = i + RANDSIZ/2 + switch i % 4 { + case 0: + a = ^(a ^ (a << 21)) + m[m2Index] + case 1: + a = (a ^ (a >> 5)) + m[m2Index] + case 2: + a = (a ^ (a << 12)) + m[m2Index] + case 3: + a = (a ^ (a >> 33)) + m[m2Index] + } + x = m[i] + y = ind(m, x) + a + b + m[i] = y + b = ind(m, y>>RANDSIZL) + x + r[i] = b + } + + // Second half + for i = RANDSIZ / 2; i < RANDSIZ; i++ { + m2Index = i - RANDSIZ/2 + switch i % 4 { + case 0: + a = ^(a ^ (a << 21)) + m[m2Index] + case 1: + a = (a ^ (a >> 5)) + m[m2Index] + case 2: + a = (a ^ (a << 12)) + m[m2Index] + case 3: + a = (a ^ (a >> 33)) + m[m2Index] + } + x = m[i] + y = ind(m, x) + a + b + m[i] = y + b = ind(m, y>>RANDSIZL) + x + r[i] = b + } + + isaac.bb = b + isaac.aa = a +} + +// Return a 64 bit random integer. +func (isaac *ISAAC) Uint64() uint64 { + if isaac.randcnt == 0 { + isaac.isaac() + isaac.randcnt = RANDSIZ + } + isaac.randcnt-- + return isaac.randrsl[isaac.randcnt] +} + +var gencycle int = 0 +var bufferFor32 uint64 = uint64(0) + +// Return a 32 bit random integer, composed of the high 32 bits of the generated 32 bit result. +func (isaac *ISAAC) Uint32() uint32 { + if gencycle == 0 { + bufferFor32 = isaac.Uint64() + gencycle = 1 + return uint32(bufferFor32 >> 32) + } + + gencycle = 0 + return uint32(bufferFor32 & 0xffffffff) +} + +// Until there is better benchmarking support in gno, you can test the performance of this PRNG with this function. +// This isn't perfect, since it will include the startup time of gno in the results, but this will give you a timing +// for generating a million random uint64 numbers on any unix based system: +// +// `time gno run -expr 'benchmarkISAAC()' isaac64.gno +func benchmarkISAAC(_iterations ...int) { + iterations := 1000000 + if len(_iterations) > 0 { + iterations = _iterations[0] + } + isaac := New() + + for i := 0; i < iterations; i++ { + _ = isaac.Uint64() + } + ufmt.Println(ufmt.Sprintf("ISAAC: generated %d uint64\n", iterations)) +} + +// The averageISAAC() function is a simple benchmarking helper to demonstrate +// the most basic statistical property of the ISAAC PRNG. +func averageISAAC(_iterations ...int) { + target := uint64(500000) + iterations := 1000000 + + ufmt.Println( + ufmt.Sprintf( + "Averaging %d random numbers. The average should be very close to %d.\n", + iterations, + target)) + + if len(_iterations) > 0 { + iterations = _iterations[0] + } + isaac := New(987654321987654321, 123456789987654321, 1, 997755331886644220) + + var average float64 = 0 + var squares []uint64 = make([]uint64, iterations) + for i := 0; i < iterations; i++ { + n := isaac.Uint64()%(target*2) + 1 + average += (float64(n) - average) / float64(i+1) + squares[i] = n + } + + sum_of_squares := uint64(0) + // transform numbers into their squares of the distance from the average + for i := 0; i < iterations; i++ { + difference := average - float64(squares[i]) + square := uint64(difference * difference) + sum_of_squares += square + } + + ufmt.Println(ufmt.Sprintf("ISAAC average of %d uint64: %f\n", iterations, average)) + ufmt.Println(ufmt.Sprintf("ISAAC standard deviation : %f\n", math.Sqrt(float64(sum_of_squares)/float64(iterations)))) + ufmt.Println(ufmt.Sprintf("ISAAC theoretical perfect deviation: %f\n", (float64(target*2)-1)/math.Sqrt(12))) +} diff --git a/examples/gno.land/p/wyhaines/rand/isaac64/isaac64_test.gno b/examples/gno.land/p/wyhaines/rand/isaac64/isaac64_test.gno new file mode 100644 index 00000000000..239e7f818fb --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/isaac64/isaac64_test.gno @@ -0,0 +1,165 @@ +package isaac64 + +import ( + "math/rand" + "testing" +) + +type OpenISAAC struct { + Randrsl [256]uint64 + Randcnt uint64 + Mm [256]uint64 + Aa, Bb, Cc uint64 + Seed [256]uint64 +} + +func TestISAACSeeding(t *testing.T) { + isaac := New() +} + +func TestISAACRand(t *testing.T) { + source := New(987654321) + rng := rand.New(source) + + // Expected outputs for the first 5 random floats with the given seed + expected := []float64{ + 0.9273376778618531, + 0.327620245173309, + 0.49315436150113456, + 0.9222536383598948, + 0.2999297342641162, + 0.4050531597269049, + 0.5321357451089953, + 0.19478000239059667, + 0.5156043950865713, + 0.9233494881511063, + } + + for i, exp := range expected { + val := rng.Float64() + if exp != val { + t.Errorf("Rand.Float64() at iteration %d: got %g, expected %g", i, val, exp) + } + } +} + +func TestISAACUint64(t *testing.T) { + isaac := New() + + expected := []uint64{ + 6781932227698873623, + 14800945299485332986, + 4114322996297394168, + 5328012296808356526, + 12789214124608876433, + 17611101631239575547, + 6877490613942924608, + 15954522518901325556, + 14180160756719376887, + 4977949063252893357, + } + + for i, exp := range expected { + val := isaac.Uint64() + if exp != val { + t.Errorf("ISAAC.Uint64() at iteration %d: got %d, expected %d", i, val, exp) + } + } +} + +func dupState(i *ISAAC) *OpenISAAC { + state := &OpenISAAC{} + state.Seed = i.seed + state.Randrsl = i.randrsl + state.Mm = i.mm + state.Aa = i.aa + state.Bb = i.bb + state.Cc = i.cc + state.Randcnt = i.randcnt + + return state +} + +func TestISAACMarshalUnmarshal(t *testing.T) { + isaac := New() + + expected1 := []uint64{ + 6781932227698873623, + 14800945299485332986, + 4114322996297394168, + 5328012296808356526, + 12789214124608876433, + } + + expected2 := []uint64{ + 17611101631239575547, + 6877490613942924608, + 15954522518901325556, + 14180160756719376887, + 4977949063252893357, + } + + for i, exp := range expected1 { + val := isaac.Uint64() + if exp != val { + t.Errorf("ISAAC.Uint64() at iteration %d: got %d, expected %d", i, val, exp) + } + } + + marshalled, err := isaac.MarshalBinary() + + t.Logf("State: [%v]\n", dupState(isaac)) + t.Logf("Marshalled State: [%x] -- %v\n", marshalled, err) + state_before := dupState(isaac) + + if err != nil { + t.Errorf("ISAAC.MarshalBinary() error: %v", err) + } + + // Advance state by one number; then check the next 5. The expectation is that they _will_ fail. + isaac.Uint64() + + for i, exp := range expected2 { + val := isaac.Uint64() + if exp == val { + t.Errorf(" Iteration %d matched %d; which is from iteration %d; something strange is happening.", (i + 6), val, (i + 5)) + } + } + + t.Logf("State before unmarshall: [%v]\n", dupState(isaac)) + + // Now restore the state of the PRNG + err = isaac.UnmarshalBinary(marshalled) + + t.Logf("State after unmarshall: [%v]\n", dupState(isaac)) + + if state_before.Seed != dupState(isaac).Seed { + t.Errorf("Seed mismatch") + } + if state_before.Randrsl != dupState(isaac).Randrsl { + t.Errorf("Randrsl mismatch") + } + if state_before.Mm != dupState(isaac).Mm { + t.Errorf("Mm mismatch") + } + if state_before.Aa != dupState(isaac).Aa { + t.Errorf("Aa mismatch") + } + if state_before.Bb != dupState(isaac).Bb { + t.Errorf("Bb mismatch") + } + if state_before.Cc != dupState(isaac).Cc { + t.Errorf("Cc mismatch") + } + if state_before.Randcnt != dupState(isaac).Randcnt { + t.Errorf("Randcnt mismatch") + } + + // Now we should be back on track for the last 5 numbers + for i, exp := range expected2 { + val := isaac.Uint64() + if exp != val { + t.Errorf("ISAAC.Uint64() at iteration %d: got %d, expected %d", (i + 5), val, exp) + } + } +} diff --git a/examples/gno.land/p/wyhaines/rand/xorshift64star/README.MD b/examples/gno.land/p/wyhaines/rand/xorshift64star/README.MD new file mode 100644 index 00000000000..00ed4412db0 --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/xorshift64star/README.MD @@ -0,0 +1,69 @@ +# package xorshift64star // import "gno.land/p/demo/math/rand/xorshift64star" + +Xorshift64* is a very fast psuedo-random number generation algorithm with strong +statistical properties. + +The default random number algorithm in gno was ported from Go's v2 rand +implementatoon, which defaults to the PCG algorithm. This algorithm is +commonly used in language PRNG implementations because it has modest seeding +requirements, and generates statistically strong randomness. + +This package provides an implementation of the Xorshift64* PRNG algorithm. +This algorithm provides strong statistical performance with most seeds (just +don't seed it with zero), and the performance of this implementation in Gno is +more than four times faster than the default PCG implementation in `math/rand`. + + +``` +Benchmark +--------- +PCG: 1000000 Uint64 generated in 15.58s +Xorshift64*: 1000000 Uint64 generated in 3.77s +Ratio: x4.11 times faster than PCG +``` + +Use it directly: + +``` +prng = xorshift64star.New() // pass a uint64 to seed it or pass nothing to seed it with entropy +``` + +Or use it as a drop-in replacement for the default PRNT in Rand: + +``` +source = xorshift64star.New() +prng := rand.New(source) +``` + +## TYPES + +``` +type Xorshift64Star struct { + // Has unexported fields. +} +``` + +Xorshift64Star is a PRNG that implements the Xorshift64* algorithm. + +`func New(seed ...uint64) *Xorshift64Star` + New() creates a new instance of the PRNG with a given seed, which should + be a uint64. If no seed is provided, the PRNG will be seeded via the + gno.land/p/demo/entropy package. + +`func (xs *Xorshift64Star) MarshalBinary() ([]byte, error)` + MarshalBinary() returns a byte array that encodes the state of the PRNG. + This can later be used with UnmarshalBinary() to restore the state of the + PRNG. MarshalBinary implements the encoding.BinaryMarshaler interface. + +`func (xs *Xorshift64Star) Seed(seed ...uint64)` + Seed() implements the rand.Source interface. It provides a way to set the + seed for the PRNG. + +`func (xs *Xorshift64Star) Uint64() uint64` + Uint64() generates the next random uint64 value. + +`func (xs *Xorshift64Star) UnmarshalBinary(data []byte) error` + UnmarshalBinary() restores the state of the PRNG from a byte array + that was created with MarshalBinary(). UnmarshalBinary implements the + encoding.BinaryUnmarshaler interface. + diff --git a/examples/gno.land/p/wyhaines/rand/xorshift64star/gno.mod b/examples/gno.land/p/wyhaines/rand/xorshift64star/gno.mod new file mode 100644 index 00000000000..bc40b1bc71b --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/xorshift64star/gno.mod @@ -0,0 +1,6 @@ +module gno.land/p/wyhaines/rand/xorshift64star + +require ( + gno.land/p/demo/entropy v0.0.0-latest + gno.land/p/demo/ufmt v0.0.0-latest +) diff --git a/examples/gno.land/p/wyhaines/rand/xorshift64star/xorshift64star.gno b/examples/gno.land/p/wyhaines/rand/xorshift64star/xorshift64star.gno new file mode 100644 index 00000000000..4934fe3a878 --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/xorshift64star/xorshift64star.gno @@ -0,0 +1,172 @@ +// Xorshift64* is a very fast psuedo-random number generation algorithm with strong +// statistical properties. +// +// The default random number algorithm in gno was ported from Go's v2 rand implementatoon, which +// defaults to the PCG algorithm. This algorithm is commonly used in language PRNG implementations +// because it has modest seeding requirements, and generates statistically strong randomness. +// +// This package provides an implementation of the Xorshift64* PRNG algorithm. This algorithm provides +// strong statistical performance with most seeds (just don't seed it with zero), and the performance +// of this implementation in Gno is more than four times faster than the default PCG implementation in +// `math/rand`. +// +// Benchmark +// --------- +// PCG: 1000000 Uint64 generated in 15.58s +// Xorshift64*: 1000000 Uint64 generated in 3.77s +// Ratio: x4.11 times faster than PCG +// +// Use it directly: +// +// prng = xorshift64star.New() // pass a uint64 to seed it or pass nothing to seed it with entropy +// +// Or use it as a drop-in replacement for the default PRNT in Rand: +// +// source = xorshift64star.New() +// prng := rand.New(source) +package xorshift64star + +import ( + "errors" + "math" + + "gno.land/p/demo/entropy" + "gno.land/p/demo/ufmt" +) + +// Xorshift64Star is a PRNG that implements the Xorshift64* algorithm. +type Xorshift64Star struct { + seed uint64 +} + +// New() creates a new instance of the PRNG with a given seed, which +// should be a uint64. If no seed is provided, the PRNG will be seeded via the +// gno.land/p/demo/entropy package. +func New(seed ...uint64) *Xorshift64Star { + xs := &Xorshift64Star{} + xs.Seed(seed...) + return xs +} + +// Seed() implements the rand.Source interface. It provides a way to set the seed for the PRNG. +func (xs *Xorshift64Star) Seed(seed ...uint64) { + if len(seed) == 0 { + e := entropy.New() + xs.seed = e.Value64() + } else { + xs.seed = seed[0] + } +} + +// beUint64() decodes a uint64 from a set of eight bytes, assuming big endian encoding. +// binary.bigEndian.Uint64, copied to avoid dependency +func beUint64(b []byte) uint64 { + _ = b[7] // bounds check hint to compiler; see golang.org/issue/14808 + return uint64(b[7]) | uint64(b[6])<<8 | uint64(b[5])<<16 | uint64(b[4])<<24 | + uint64(b[3])<<32 | uint64(b[2])<<40 | uint64(b[1])<<48 | uint64(b[0])<<56 +} + +// bePutUint64() encodes a uint64 into a buffer of eight bytes. +// binary.bigEndian.PutUint64, copied to avoid dependency +func bePutUint64(b []byte, v uint64) { + _ = b[7] // early bounds check to guarantee safety of writes below + b[0] = byte(v >> 56) + b[1] = byte(v >> 48) + b[2] = byte(v >> 40) + b[3] = byte(v >> 32) + b[4] = byte(v >> 24) + b[5] = byte(v >> 16) + b[6] = byte(v >> 8) + b[7] = byte(v) +} + +// A label to identify the marshalled data. +var marshalXorshift64StarLabel = []byte("xorshift64*:") + +// MarshalBinary() returns a byte array that encodes the state of the PRNG. This can later be used +// with UnmarshalBinary() to restore the state of the PRNG. +// MarshalBinary implements the encoding.BinaryMarshaler interface. +func (xs *Xorshift64Star) MarshalBinary() ([]byte, error) { + b := make([]byte, 20) + copy(b, marshalXorshift64StarLabel) + bePutUint64(b[12:], xs.seed) + return b, nil +} + +// errUnmarshalXorshift64Star is returned when unmarshalling fails. +var errUnmarshalXorshift64Star = errors.New("invalid Xorshift64* encoding") + +// UnmarshalBinary() restores the state of the PRNG from a byte array that was created with MarshalBinary(). +// UnmarshalBinary implements the encoding.BinaryUnmarshaler interface. +func (xs *Xorshift64Star) UnmarshalBinary(data []byte) error { + if len(data) != 20 || string(data[:12]) != string(marshalXorshift64StarLabel) { + return errUnmarshalXorshift64Star + } + xs.seed = beUint64(data[12:]) + return nil +} + +// Uint64() generates the next random uint64 value. +func (xs *Xorshift64Star) Uint64() uint64 { + xs.seed ^= xs.seed >> 12 + xs.seed ^= xs.seed << 25 + xs.seed ^= xs.seed >> 27 + xs.seed *= 2685821657736338717 + return xs.seed // Operations naturally wrap around in uint64 +} + +// Until there is better benchmarking support in gno, you can test the performance of this PRNG with this function. +// This isn't perfect, since it will include the startup time of gno in the results, but this will give you a timing +// for generating a million random uint64 numbers on any unix based system: +// +// `time gno run -expr 'benchmarkXorshift64Star()' xorshift64star.gno +func benchmarkXorshift64Star(_iterations ...int) { + iterations := 1000000 + if len(_iterations) > 0 { + iterations = _iterations[0] + } + xs64s := New() + + for i := 0; i < iterations; i++ { + _ = xs64s.Uint64() + } + ufmt.Println(ufmt.Sprintf("Xorshift64*: generate %d uint64\n", iterations)) +} + +// The averageXorshift64Star() function is a simple benchmarking helper to demonstrate +// the most basic statistical property of the Xorshift64* PRNG. +func averageXorshift64Star(_iterations ...int) { + target := uint64(500000) + iterations := 1000000 + var squares [1000000]uint64 + + ufmt.Println( + ufmt.Sprintf( + "Averaging %d random numbers. The average should be very close to %d.\n", + iterations, + target)) + + if len(_iterations) > 0 { + iterations = _iterations[0] + } + xs64s := New() + + var average float64 = 0 + for i := 0; i < iterations; i++ { + n := xs64s.Uint64()%(target*2) + 1 + average += (float64(n) - average) / float64(i+1) + squares[i] = n + } + + sum_of_squares := uint64(0) + // transform numbers into their squares of the distance from the average + for i := 0; i < iterations; i++ { + difference := average - float64(squares[i]) + square := uint64(difference * difference) + sum_of_squares += square + } + + ufmt.Println(ufmt.Sprintf("Xorshift64* average of %d uint64: %f\n", iterations, average)) + ufmt.Println(ufmt.Sprintf("Xorshift64* standard deviation : %f\n", math.Sqrt(float64(sum_of_squares)/float64(iterations)))) + ufmt.Println(ufmt.Sprintf("Xorshift64* theoretical perfect deviation: %f\n", (float64(target*2)-1)/math.Sqrt(12))) +} diff --git a/examples/gno.land/p/wyhaines/rand/xorshift64star/xorshift64star_test.gno b/examples/gno.land/p/wyhaines/rand/xorshift64star/xorshift64star_test.gno new file mode 100644 index 00000000000..8a73bd9718d --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/xorshift64star/xorshift64star_test.gno @@ -0,0 +1,134 @@ +package xorshift64star + +import ( + "math/rand" + "testing" +) + +func TestXorshift64StarSeeding(t *testing.T) { + xs64s := New() + value1 := xs64s.Uint64() + + xs64s = New(987654321) + value2 := xs64s.Uint64() + + if value1 != 5083824587905981259 || value2 != 18211065302896784785 || value1 == value2 { + t.Errorf("Expected 5083824587905981259 to be != to 18211065302896784785; got: %d == %d", value1, value2) + } +} + +func TestXorshift64StarRand(t *testing.T) { + source := New(987654321) + rng := rand.New(source) + + // Expected outputs for the first 5 random floats with the given seed + expected := []float64{ + .8344002228310946, + 0.01777174153236205, + 0.23521769507865276, + 0.5387610198576143, + 0.631539862225968, + 0.9369068148346704, + 0.6387002315083188, + 0.5047507613688854, + 0.5208486273732391, + 0.25023746271541747, + } + + for i, exp := range expected { + val := rng.Float64() + if exp != val { + t.Errorf("Rand.Float64() at iteration %d: got %g, expected %g", i, val, exp) + } + } +} + +func TestXorshift64StarUint64(t *testing.T) { + xs64s := New() + + expected := []uint64{ + 5083824587905981259, + 4607286371009545754, + 2070557085263023674, + 14094662988579565368, + 2910745910478213381, + 18037409026311016155, + 17169624916429864153, + 10459214929523155306, + 11840179828060641081, + 1198750959721587199, + } + + for i, exp := range expected { + val := xs64s.Uint64() + if exp != val { + t.Errorf("Xorshift64Star.Uint64() at iteration %d: got %d, expected %d", i, val, exp) + } + } +} + +func TestXorshift64StarMarshalUnmarshal(t *testing.T) { + xs64s := New() + + expected1 := []uint64{ + 5083824587905981259, + 4607286371009545754, + 2070557085263023674, + 14094662988579565368, + 2910745910478213381, + } + + expected2 := []uint64{ + 18037409026311016155, + 17169624916429864153, + 10459214929523155306, + 11840179828060641081, + 1198750959721587199, + } + + for i, exp := range expected1 { + val := xs64s.Uint64() + if exp != val { + t.Errorf("Xorshift64Star.Uint64() at iteration %d: got %d, expected %d", i, val, exp) + } + } + + marshalled, err := xs64s.MarshalBinary() + + t.Logf("Original State: [%x]\n", xs64s.seed) + t.Logf("Marshalled State: [%x] -- %v\n", marshalled, err) + state_before := xs64s.seed + + if err != nil { + t.Errorf("Xorshift64Star.MarshalBinary() error: %v", err) + } + + // Advance state by one number; then check the next 5. The expectation is that they _will_ fail. + xs64s.Uint64() + + for i, exp := range expected2 { + val := xs64s.Uint64() + if exp == val { + t.Errorf(" Iteration %d matched %d; which is from iteration %d; something strange is happening.", (i + 6), val, (i + 5)) + } + } + + t.Logf("State before unmarshall: [%x]\n", xs64s.seed) + + // Now restore the state of the PRNG + err = xs64s.UnmarshalBinary(marshalled) + + t.Logf("State after unmarshall: [%x]\n", xs64s.seed) + + if state_before != xs64s.seed { + t.Errorf("States before and after marshal/unmarshal are not equal; go %x and %x", state_before, xs64s.seed) + } + + // Now we should be back on track for the last 5 numbers + for i, exp := range expected2 { + val := xs64s.Uint64() + if exp != val { + t.Errorf("Xorshift64Star.Uint64() at iteration %d: got %d, expected %d", (i + 5), val, exp) + } + } +} diff --git a/examples/gno.land/p/wyhaines/rand/xorshiftr128plus/README.MD b/examples/gno.land/p/wyhaines/rand/xorshiftr128plus/README.MD new file mode 100644 index 00000000000..444d1e1cdd9 --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/xorshiftr128plus/README.MD @@ -0,0 +1,60 @@ +# package xorshiftr128plus // import "gno.land/p/demo/math/rand/xorshiftr128plus" + +Xorshiftr128+ is a very fast psuedo-random number generation algorithm with +strong statistical properties. + +The default random number algorithm in gno was ported from Go's v2 rand +implementatoon, which defaults to the PCG algorithm. This algorithm is +commonly used in language PRNG implementations because it has modest seeding +requirements, and generates statistically strong randomness. + +This package provides an implementation of the Xorshiftr128+ PRNG algorithm. +This algorithm provides strong statistical performance with most seeds (just +don't seed it with zeros), and the performance of this implementation in Gno is +more than four times faster than the default PCG implementation in `math/rand`. + +``` +Benchmark +--------- +PCG: 1000000 Uint64 generated in 15.48s +Xorshiftr128+: 1000000 Uint64 generated in 3.22s +Ratio: x4.81 times faster than PCG +``` + +Use it directly: + +``` +prng = xorshiftr128plus.New() // pass a uint64 to seed it or pass nothing to seed it with entropy +``` + +Or use it as a drop-in replacement for the default PRNT in Rand: + +``` +source = xorshiftr128plus.New() +prng := rand.New(source) +``` + +## TYPES + +``` +type Xorshiftr128Plus struct { + // Has unexported fields. +} +``` + +`func New(seeds ...uint64) *Xorshiftr128Plus` + +`func (xs *Xorshiftr128Plus) MarshalBinary() ([]byte, error)` + MarshalBinary() returns a byte array that encodes the state of the PRNG. + This can later be used with UnmarshalBinary() to restore the state of the + PRNG. MarshalBinary implements the encoding.BinaryMarshaler interface. + +`func (x *Xorshiftr128Plus) Seed(s1, s2 uint64)` + +`func (x *Xorshiftr128Plus) Uint64() uint64` + +`func (xs *Xorshiftr128Plus) UnmarshalBinary(data []byte) error` + UnmarshalBinary() restores the state of the PRNG from a byte array + that was created with MarshalBinary(). UnmarshalBinary implements the + encoding.BinaryUnmarshaler interface. + diff --git a/examples/gno.land/p/wyhaines/rand/xorshiftr128plus/gno.mod b/examples/gno.land/p/wyhaines/rand/xorshiftr128plus/gno.mod new file mode 100644 index 00000000000..c778fc72550 --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/xorshiftr128plus/gno.mod @@ -0,0 +1,6 @@ +module gno.land/p/wyhaines/rand/xorshiftr128plus + +require ( + gno.land/p/demo/entropy v0.0.0-latest + gno.land/p/demo/ufmt v0.0.0-latest +) diff --git a/examples/gno.land/p/wyhaines/rand/xorshiftr128plus/xorshiftr128plus.gno b/examples/gno.land/p/wyhaines/rand/xorshiftr128plus/xorshiftr128plus.gno new file mode 100644 index 00000000000..d950ab5108a --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/xorshiftr128plus/xorshiftr128plus.gno @@ -0,0 +1,186 @@ +// Xorshiftr128+ is a very fast psuedo-random number generation algorithm with strong +// statistical properties. +// +// The default random number algorithm in gno was ported from Go's v2 rand implementatoon, which +// defaults to the PCG algorithm. This algorithm is commonly used in language PRNG implementations +// because it has modest seeding requirements, and generates statistically strong randomness. +// +// This package provides an implementation of the Xorshiftr128+ PRNG algorithm. This algorithm provides +// strong statistical performance with most seeds (just don't seed it with zeros), and the performance +// of this implementation in Gno is more than four times faster than the default PCG implementation in +// `math/rand`. +// +// Benchmark +// --------- +// PCG: 1000000 Uint64 generated in 15.48s +// Xorshiftr128+: 1000000 Uint64 generated in 3.22s +// Ratio: x4.81 times faster than PCG +// +// Use it directly: +// +// prng = xorshiftr128plus.New() // pass a uint64 to seed it or pass nothing to seed it with entropy +// +// Or use it as a drop-in replacement for the default PRNT in Rand: +// +// source = xorshiftr128plus.New() +// prng := rand.New(source) +package xorshiftr128plus + +import ( + "errors" + "math" + + "gno.land/p/demo/entropy" + "gno.land/p/demo/ufmt" +) + +type Xorshiftr128Plus struct { + seed [2]uint64 // Seeds +} + +func New(seeds ...uint64) *Xorshiftr128Plus { + var s1, s2 uint64 + seed_length := len(seeds) + if seed_length < 2 { + e := entropy.New() + if seed_length == 0 { + s1 = e.Value64() + s2 = e.Value64() + } else { + s1 = seeds[0] + s2 = e.Value64() + } + } else { + s1 = seeds[0] + s2 = seeds[1] + } + + prng := &Xorshiftr128Plus{} + prng.Seed(s1, s2) + return prng +} + +func (x *Xorshiftr128Plus) Seed(s1, s2 uint64) { + if s1 == 0 && s2 == 0 { + panic("Seeds must not both be zero") + } + x.seed[0] = s1 + x.seed[1] = s2 +} + +// beUint64() decodes a uint64 from a set of eight bytes, assuming big endian encoding. +// binary.bigEndian.Uint64, copied to avoid dependency +func beUint64(b []byte) uint64 { + _ = b[7] // bounds check hint to compiler; see golang.org/issue/14808 + return uint64(b[7]) | uint64(b[6])<<8 | uint64(b[5])<<16 | uint64(b[4])<<24 | + uint64(b[3])<<32 | uint64(b[2])<<40 | uint64(b[1])<<48 | uint64(b[0])<<56 +} + +// bePutUint64() encodes a uint64 into a buffer of eight bytes. +// binary.bigEndian.PutUint64, copied to avoid dependency +func bePutUint64(b []byte, v uint64) { + _ = b[7] // early bounds check to guarantee safety of writes below + b[0] = byte(v >> 56) + b[1] = byte(v >> 48) + b[2] = byte(v >> 40) + b[3] = byte(v >> 32) + b[4] = byte(v >> 24) + b[5] = byte(v >> 16) + b[6] = byte(v >> 8) + b[7] = byte(v) +} + +// A label to identify the marshalled data. +var marshalXorshiftr128PlusLabel = []byte("xorshiftr128+:") + +// MarshalBinary() returns a byte array that encodes the state of the PRNG. This can later be used +// with UnmarshalBinary() to restore the state of the PRNG. +// MarshalBinary implements the encoding.BinaryMarshaler interface. +func (xs *Xorshiftr128Plus) MarshalBinary() ([]byte, error) { + b := make([]byte, 30) + copy(b, marshalXorshiftr128PlusLabel) + bePutUint64(b[14:], xs.seed[0]) + bePutUint64(b[22:], xs.seed[1]) + return b, nil +} + +// errUnmarshalXorshiftr128Plus is returned when unmarshalling fails. +var errUnmarshalXorshiftr128Plus = errors.New("invalid Xorshiftr128Plus encoding") + +// UnmarshalBinary() restores the state of the PRNG from a byte array that was created with MarshalBinary(). +// UnmarshalBinary implements the encoding.BinaryUnmarshaler interface. +func (xs *Xorshiftr128Plus) UnmarshalBinary(data []byte) error { + if len(data) != 30 || string(data[:14]) != string(marshalXorshiftr128PlusLabel) { + return errUnmarshalXorshiftr128Plus + } + xs.seed[0] = beUint64(data[14:]) + xs.seed[1] = beUint64(data[22:]) + return nil +} + +func (x *Xorshiftr128Plus) Uint64() uint64 { + x0 := x.seed[0] + x1 := x.seed[1] + x.seed[0] = x1 + x0 ^= x0 << 23 + x0 ^= x0 >> 17 + x0 ^= x1 + x.seed[1] = x0 + x1 + return x.seed[1] +} + +// Until there is better benchmarking support in gno, you can test the performance of this PRNG with this function. +// This isn't perfect, since it will include the startup time of gno in the results, but this will give you a timing +// for generating a million random uint64 numbers on any unix based system: +// +// `time gno run -expr 'benchmarkXorshiftr128Plus()' xorshiftr128plus.gno +func benchmarkXorshiftr128Plus(_iterations ...int) { + iterations := 1000000 + if len(_iterations) > 0 { + iterations = _iterations[0] + } + xs128p := New() + + for i := 0; i < iterations; i++ { + _ = xs128p.Uint64() + } + ufmt.Println(ufmt.Sprintf("Xorshiftr128Plus: generate %d uint64\n", iterations)) +} + +// The averageXorshiftr128Plus() function is a simple benchmarking helper to demonstrate +// the most basic statistical property of the Xorshiftr128+ PRNG. +func averageXorshiftr128Plus(_iterations ...int) { + target := uint64(500000) + iterations := 1000000 + var squares [1000000]uint64 + + ufmt.Println( + ufmt.Sprintf( + "Averaging %d random numbers. The average should be very close to %d.\n", + iterations, + target)) + + if len(_iterations) > 0 { + iterations = _iterations[0] + } + xs128p := New() + + var average float64 = 0 + for i := 0; i < iterations; i++ { + n := xs128p.Uint64()%(target*2) + 1 + average += (float64(n) - average) / float64(i+1) + squares[i] = n + } + + sum_of_squares := uint64(0) + // transform numbers into their squares of the distance from the average + for i := 0; i < iterations; i++ { + difference := average - float64(squares[i]) + square := uint64(difference * difference) + sum_of_squares += square + } + + ufmt.Println(ufmt.Sprintf("Xorshiftr128+ average of %d uint64: %f\n", iterations, average)) + ufmt.Println(ufmt.Sprintf("Xorshiftr128+ standard deviation : %f\n", math.Sqrt(float64(sum_of_squares)/float64(iterations)))) + ufmt.Println(ufmt.Sprintf("Xorshiftr128+ theoretical perfect deviation: %f\n", (float64(target*2)-1)/math.Sqrt(12))) +} diff --git a/examples/gno.land/p/wyhaines/rand/xorshiftr128plus/xorshiftr128plus_test.gno b/examples/gno.land/p/wyhaines/rand/xorshiftr128plus/xorshiftr128plus_test.gno new file mode 100644 index 00000000000..c5d86edd073 --- /dev/null +++ b/examples/gno.land/p/wyhaines/rand/xorshiftr128plus/xorshiftr128plus_test.gno @@ -0,0 +1,142 @@ +package xorshiftr128plus + +import ( + "math/rand" + "testing" +) + +func TestXorshift64StarSeeding(t *testing.T) { + xs128p := New() + value1 := xs128p.Uint64() + + xs128p = New(987654321) + value2 := xs128p.Uint64() + + xs128p = New(987654321, 9876543210) + value3 := xs128p.Uint64() + + if value1 != 13970141264473760763 || + value2 != 17031892808144362974 || + value3 != 8285073084540510 || + value1 == value2 || + value2 == value3 || + value1 == value3 { + t.Errorf("Expected three different values: 13970141264473760763, 17031892808144362974, and 8285073084540510\n got: %d, %d, %d", value1, value2, value3) + } +} + +func TestXorshiftr128PlusRand(t *testing.T) { + source := New(987654321) + rng := rand.New(source) + + // Expected outputs for the first 5 random floats with the given seed + expected := []float64{ + 0.9199548549485674, + 0.0027491282372705816, + 0.31493362274701164, + 0.3531250819119609, + 0.09957852858060356, + 0.731941362705936, + 0.3476937688876708, + 0.1444018086140385, + 0.9106467321832331, + 0.8024870151488901, + } + + for i, exp := range expected { + val := rng.Float64() + if exp != val { + t.Errorf("Rand.Float64() at iteration %d: got %g, expected %g", i, val, exp) + } + } +} + +func TestXorshiftr128PlusUint64(t *testing.T) { + xs128p := New(987654321, 9876543210) + + expected := []uint64{ + 8285073084540510, + 97010855169053386, + 11353359435625603792, + 10289232744262291728, + 14019961444418950453, + 15829492476941720545, + 2764732928842099222, + 6871047144273883379, + 16142204260470661970, + 11803223757041229095, + } + + for i, exp := range expected { + val := xs128p.Uint64() + if exp != val { + t.Errorf("Xorshiftr128Plus.Uint64() at iteration %d: got %d, expected %d", i, val, exp) + } + } +} + +func TestXorshiftr128PlusMarshalUnmarshal(t *testing.T) { + xs128p := New(987654321, 9876543210) + + expected1 := []uint64{ + 8285073084540510, + 97010855169053386, + 11353359435625603792, + 10289232744262291728, + 14019961444418950453, + } + + expected2 := []uint64{ + 15829492476941720545, + 2764732928842099222, + 6871047144273883379, + 16142204260470661970, + 11803223757041229095, + } + + for i, exp := range expected1 { + val := xs128p.Uint64() + if exp != val { + t.Errorf("Xorshiftr128Plus.Uint64() at iteration %d: got %d, expected %d", i, val, exp) + } + } + + marshalled, err := xs128p.MarshalBinary() + + t.Logf("Original State: [%x]\n", xs128p.seed) + t.Logf("Marshalled State: [%x] -- %v\n", marshalled, err) + state_before := xs128p.seed + + if err != nil { + t.Errorf("Xorshiftr128Plus.MarshalBinary() error: %v", err) + } + + // Advance state by one number; then check the next 5. The expectation is that they _will_ fail. + xs128p.Uint64() + + for i, exp := range expected2 { + val := xs128p.Uint64() + if exp == val { + t.Errorf(" Iteration %d matched %d; which is from iteration %d; something strange is happening.", (i + 6), val, (i + 5)) + } + } + + t.Logf("State before unmarshall: [%x]\n", xs128p.seed) + + // Now restore the state of the PRNG + err = xs128p.UnmarshalBinary(marshalled) + + t.Logf("State after unmarshall: [%x]\n", xs128p.seed) + + if state_before != xs128p.seed { + t.Errorf("States before and after marshal/unmarshal are not equal; go %x and %x", state_before, xs128p.seed) + } + + // Now we should be back on track for the last 5 numbers + for i, exp := range expected2 { + val := xs128p.Uint64() + if exp != val { + t.Errorf("Xorshiftr128Plus.Uint64() at iteration %d: got %d, expected %d", (i + 5), val, exp) + } + } +}