xorshift64star.gno
5.83 Kb ยท 172 lines
1// Xorshift64* 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 Xorshift64* PRNG algorithm. This algorithm provides
9// strong statistical performance with most seeds (just don't seed it with zero), 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.58s
16// Xorshift64*: 1000000 Uint64 generated in 3.77s
17// Ratio: x4.11 times faster than PCG
18//
19// Use it directly:
20//
21// prng = xorshift64star.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 = xorshift64star.New()
26// prng := rand.New(source)
27package xorshift64star
28
29import (
30 "errors"
31 "math"
32
33 "gno.land/p/demo/entropy"
34 "gno.land/p/nt/ufmt"
35)
36
37// Xorshift64Star is a PRNG that implements the Xorshift64* algorithm.
38type Xorshift64Star struct {
39 seed uint64
40}
41
42// New() creates a new instance of the PRNG with a given seed, which
43// should be a uint64. If no seed is provided, the PRNG will be seeded via the
44// gno.land/p/demo/entropy package.
45func New(seed ...uint64) *Xorshift64Star {
46 xs := &Xorshift64Star{}
47 xs.Seed(seed...)
48 return xs
49}
50
51// Seed() implements the rand.Source interface. It provides a way to set the seed for the PRNG.
52func (xs *Xorshift64Star) Seed(seed ...uint64) {
53 if len(seed) == 0 {
54 e := entropy.New()
55 xs.seed = e.Value64()
56 } else {
57 xs.seed = seed[0]
58 }
59}
60
61// beUint64() decodes a uint64 from a set of eight bytes, assuming big endian encoding.
62// binary.bigEndian.Uint64, copied to avoid dependency
63func beUint64(b []byte) uint64 {
64 _ = b[7] // bounds check hint to compiler; see golang.org/issue/14808
65 return uint64(b[7]) | uint64(b[6])<<8 | uint64(b[5])<<16 | uint64(b[4])<<24 |
66 uint64(b[3])<<32 | uint64(b[2])<<40 | uint64(b[1])<<48 | uint64(b[0])<<56
67}
68
69// bePutUint64() encodes a uint64 into a buffer of eight bytes.
70// binary.bigEndian.PutUint64, copied to avoid dependency
71func bePutUint64(b []byte, v uint64) {
72 _ = b[7] // early bounds check to guarantee safety of writes below
73 b[0] = byte(v >> 56)
74 b[1] = byte(v >> 48)
75 b[2] = byte(v >> 40)
76 b[3] = byte(v >> 32)
77 b[4] = byte(v >> 24)
78 b[5] = byte(v >> 16)
79 b[6] = byte(v >> 8)
80 b[7] = byte(v)
81}
82
83// A label to identify the marshalled data.
84var marshalXorshift64StarLabel = []byte("xorshift64*:")
85
86// MarshalBinary() returns a byte array that encodes the state of the PRNG. This can later be used
87// with UnmarshalBinary() to restore the state of the PRNG.
88// MarshalBinary implements the encoding.BinaryMarshaler interface.
89func (xs *Xorshift64Star) MarshalBinary() ([]byte, error) {
90 b := make([]byte, 20)
91 copy(b, marshalXorshift64StarLabel)
92 bePutUint64(b[12:], xs.seed)
93 return b, nil
94}
95
96// errUnmarshalXorshift64Star is returned when unmarshalling fails.
97var errUnmarshalXorshift64Star = errors.New("invalid Xorshift64* encoding")
98
99// UnmarshalBinary() restores the state of the PRNG from a byte array that was created with MarshalBinary().
100// UnmarshalBinary implements the encoding.BinaryUnmarshaler interface.
101func (xs *Xorshift64Star) UnmarshalBinary(data []byte) error {
102 if len(data) != 20 || string(data[:12]) != string(marshalXorshift64StarLabel) {
103 return errUnmarshalXorshift64Star
104 }
105 xs.seed = beUint64(data[12:])
106 return nil
107}
108
109// Uint64() generates the next random uint64 value.
110func (xs *Xorshift64Star) Uint64() uint64 {
111 xs.seed ^= xs.seed >> 12
112 xs.seed ^= xs.seed << 25
113 xs.seed ^= xs.seed >> 27
114 xs.seed *= 2685821657736338717
115 return xs.seed // Operations naturally wrap around in uint64
116}
117
118// Until there is better benchmarking support in gno, you can test the performance of this PRNG with this function.
119// This isn't perfect, since it will include the startup time of gno in the results, but this will give you a timing
120// for generating a million random uint64 numbers on any unix based system:
121//
122// `time gno run -expr 'benchmarkXorshift64Star()' xorshift64star.gno
123func benchmarkXorshift64Star(_iterations ...int) {
124 iterations := 1000000
125 if len(_iterations) > 0 {
126 iterations = _iterations[0]
127 }
128 xs64s := New()
129
130 for i := 0; i < iterations; i++ {
131 _ = xs64s.Uint64()
132 }
133 ufmt.Println(ufmt.Sprintf("Xorshift64*: generate %d uint64\n", iterations))
134}
135
136// The averageXorshift64Star() function is a simple benchmarking helper to demonstrate
137// the most basic statistical property of the Xorshift64* PRNG.
138func averageXorshift64Star(_iterations ...int) {
139 target := uint64(500000)
140 iterations := 1000000
141 var squares [1000000]uint64
142
143 ufmt.Println(
144 ufmt.Sprintf(
145 "Averaging %d random numbers. The average should be very close to %d.\n",
146 iterations,
147 target))
148
149 if len(_iterations) > 0 {
150 iterations = _iterations[0]
151 }
152 xs64s := New()
153
154 var average float64 = 0
155 for i := 0; i < iterations; i++ {
156 n := xs64s.Uint64()%(target*2) + 1
157 average += (float64(n) - average) / float64(i+1)
158 squares[i] = n
159 }
160
161 sum_of_squares := uint64(0)
162 // transform numbers into their squares of the distance from the average
163 for i := 0; i < iterations; i++ {
164 difference := average - float64(squares[i])
165 square := uint64(difference * difference)
166 sum_of_squares += square
167 }
168
169 ufmt.Println(ufmt.Sprintf("Xorshift64* average of %d uint64: %f\n", iterations, average))
170 ufmt.Println(ufmt.Sprintf("Xorshift64* standard deviation : %f\n", math.Sqrt(float64(sum_of_squares)/float64(iterations))))
171 ufmt.Println(ufmt.Sprintf("Xorshift64* theoretical perfect deviation: %f\n", (float64(target*2)-1)/math.Sqrt(12)))
172}
