xorshiftr128plus.gno

  1// Xorshiftr128+ is a very fast psuedo-random number generation algorithm with strong
  2// statistical properties.
  3//
  4// The default random number algorithm in gno was ported from Go's v2 rand implementatoon, which
  5// defaults to the PCG algorithm. This algorithm is commonly used in language PRNG implementations
  6// because it has modest seeding requirements, and generates statistically strong randomness.
  7//
  8// This package provides an implementation of the Xorshiftr128+ PRNG algorithm. This algorithm provides
  9// strong statistical performance with most seeds (just don't seed it with zeros), and the performance
 10// of this implementation in Gno is more than four times faster than the default PCG implementation in
 11// `math/rand`.
 12//
 13//	Benchmark
 14//	---------
 15//	PCG:           1000000 Uint64 generated in 15.48s
 16//	Xorshiftr128+: 1000000 Uint64 generated in 3.22s
 17//	Ratio:         x4.81 times faster than PCG
 18//
 19// Use it directly:
 20//
 21//	prng = xorshiftr128plus.New() // pass a uint64 to seed it or pass nothing to seed it with entropy
 22//
 23// Or use it as a drop-in replacement for the default PRNT in Rand:
 24//
 25//	source = xorshiftr128plus.New()
 26//	prng := rand.New(source)
 27package xorshiftr128plus
 28
 29import (
 30	"errors"
 31	"math"
 32
 33	"gno.land/p/demo/entropy"
 34	"gno.land/p/nt/ufmt"
 35)
 36
 37type Xorshiftr128Plus struct {
 38	seed [2]uint64 // Seeds
 39}
 40
 41func New(seeds ...uint64) *Xorshiftr128Plus {
 42	var s1, s2 uint64
 43	seed_length := len(seeds)
 44	if seed_length < 2 {
 45		e := entropy.New()
 46		if seed_length == 0 {
 47			s1 = e.Value64()
 48			s2 = e.Value64()
 49		} else {
 50			s1 = seeds[0]
 51			s2 = e.Value64()
 52		}
 53	} else {
 54		s1 = seeds[0]
 55		s2 = seeds[1]
 56	}
 57
 58	prng := &Xorshiftr128Plus{}
 59	prng.Seed(s1, s2)
 60	return prng
 61}
 62
 63func (x *Xorshiftr128Plus) Seed(s1, s2 uint64) {
 64	if s1 == 0 && s2 == 0 {
 65		panic("Seeds must not both be zero")
 66	}
 67	x.seed[0] = s1
 68	x.seed[1] = s2
 69}
 70
 71// beUint64() decodes a uint64 from a set of eight bytes, assuming big endian encoding.
 72// binary.bigEndian.Uint64, copied to avoid dependency
 73func beUint64(b []byte) uint64 {
 74	_ = b[7] // bounds check hint to compiler; see golang.org/issue/14808
 75	return uint64(b[7]) | uint64(b[6])<<8 | uint64(b[5])<<16 | uint64(b[4])<<24 |
 76		uint64(b[3])<<32 | uint64(b[2])<<40 | uint64(b[1])<<48 | uint64(b[0])<<56
 77}
 78
 79// bePutUint64() encodes a uint64 into a buffer of eight bytes.
 80// binary.bigEndian.PutUint64, copied to avoid dependency
 81func bePutUint64(b []byte, v uint64) {
 82	_ = b[7] // early bounds check to guarantee safety of writes below
 83	b[0] = byte(v >> 56)
 84	b[1] = byte(v >> 48)
 85	b[2] = byte(v >> 40)
 86	b[3] = byte(v >> 32)
 87	b[4] = byte(v >> 24)
 88	b[5] = byte(v >> 16)
 89	b[6] = byte(v >> 8)
 90	b[7] = byte(v)
 91}
 92
 93// A label to identify the marshalled data.
 94var marshalXorshiftr128PlusLabel = []byte("xorshiftr128+:")
 95
 96// MarshalBinary() returns a byte array that encodes the state of the PRNG. This can later be used
 97// with UnmarshalBinary() to restore the state of the PRNG.
 98// MarshalBinary implements the encoding.BinaryMarshaler interface.
 99func (xs *Xorshiftr128Plus) MarshalBinary() ([]byte, error) {
100	b := make([]byte, 30)
101	copy(b, marshalXorshiftr128PlusLabel)
102	bePutUint64(b[14:], xs.seed[0])
103	bePutUint64(b[22:], xs.seed[1])
104	return b, nil
105}
106
107// errUnmarshalXorshiftr128Plus is returned when unmarshalling fails.
108var errUnmarshalXorshiftr128Plus = errors.New("invalid Xorshiftr128Plus encoding")
109
110// UnmarshalBinary() restores the state of the PRNG from a byte array that was created with MarshalBinary().
111// UnmarshalBinary implements the encoding.BinaryUnmarshaler interface.
112func (xs *Xorshiftr128Plus) UnmarshalBinary(data []byte) error {
113	if len(data) != 30 || string(data[:14]) != string(marshalXorshiftr128PlusLabel) {
114		return errUnmarshalXorshiftr128Plus
115	}
116	xs.seed[0] = beUint64(data[14:])
117	xs.seed[1] = beUint64(data[22:])
118	return nil
119}
120
121func (x *Xorshiftr128Plus) Uint64() uint64 {
122	x0 := x.seed[0]
123	x1 := x.seed[1]
124	x.seed[0] = x1
125	x0 ^= x0 << 23
126	x0 ^= x0 >> 17
127	x0 ^= x1
128	x.seed[1] = x0 + x1
129	return x.seed[1]
130}
131
132// Until there is better benchmarking support in gno, you can test the performance of this PRNG with this function.
133// This isn't perfect, since it will include the startup time of gno in the results, but this will give you a timing
134// for generating a million random uint64 numbers on any unix based system:
135//
136// `time gno run -expr 'benchmarkXorshiftr128Plus()' xorshiftr128plus.gno
137func benchmarkXorshiftr128Plus(_iterations ...int) {
138	iterations := 1000000
139	if len(_iterations) > 0 {
140		iterations = _iterations[0]
141	}
142	xs128p := New()
143
144	for i := 0; i < iterations; i++ {
145		_ = xs128p.Uint64()
146	}
147	ufmt.Println(ufmt.Sprintf("Xorshiftr128Plus: generate %d uint64\n", iterations))
148}
149
150// The averageXorshiftr128Plus() function is a simple benchmarking helper to demonstrate
151// the most basic statistical property of the Xorshiftr128+ PRNG.
152func averageXorshiftr128Plus(_iterations ...int) {
153	target := uint64(500000)
154	iterations := 1000000
155	var squares [1000000]uint64
156
157	ufmt.Println(
158		ufmt.Sprintf(
159			"Averaging %d random numbers. The average should be very close to %d.\n",
160			iterations,
161			target))
162
163	if len(_iterations) > 0 {
164		iterations = _iterations[0]
165	}
166	xs128p := New()
167
168	var average float64 = 0
169	for i := 0; i < iterations; i++ {
170		n := xs128p.Uint64()%(target*2) + 1
171		average += (float64(n) - average) / float64(i+1)
172		squares[i] = n
173	}
174
175	sum_of_squares := uint64(0)
176	// transform numbers into their squares of the distance from the average
177	for i := 0; i < iterations; i++ {
178		difference := average - float64(squares[i])
179		square := uint64(difference * difference)
180		sum_of_squares += square
181	}
182
183	ufmt.Println(ufmt.Sprintf("Xorshiftr128+ average of %d uint64: %f\n", iterations, average))
184	ufmt.Println(ufmt.Sprintf("Xorshiftr128+ standard deviation  : %f\n", math.Sqrt(float64(sum_of_squares)/float64(iterations))))
185	ufmt.Println(ufmt.Sprintf("Xorshiftr128+ theoretical perfect deviation: %f\n", (float64(target*2)-1)/math.Sqrt(12)))
186}