isaac64.gno
11.12 Kb ยท 429 lines
1// This is a port of the 64-bit version of the ISAAC cryptographically secure PRNG, originally
2// based on the reference implementation found at https://burtleburtle.net/bob/rand/isaacafa.html
3//
4// ISAAC has excellent statistical properties, with long cycle times, and uniformly distributed,
5// unbiased, and unpredictable number generation. It can not be distinguished from real random
6// data, and in three decades of scrutiny, no practical attacks have been found.
7//
8// The default random number algorithm in gno was ported from Go's v2 rand implementatoon, which
9// defaults to the PCG algorithm. This algorithm is commonly used in language PRNG implementations
10// because it has modest seeding requirements, and generates statistically strong randomness.
11//
12// This package provides an implementation of the 64-bit ISAAC PRNG algorithm. This algorithm
13// provides very strong statistical performance, and is cryptographically secure, while still
14// being substantially faster than the default PCG implementation in `math/rand`.
15//
16// Note that the approach to seeing with ISAAC is very important for best results, and seeding with
17// ISAAC is not as simple as seeding with a single uint64 value. The ISAAC algorithm requires a
18// 256-element seed. If used for cryptographic purposes, this will likely require entropy generated
19// off-chain for actual cryptographically secure seeding. For other purposes, however, one can
20// utilize the built-in seeding mechanism, which will leverage the xorshiftr128plus PRNG to generate
21// any missing seeds if fewer than 256 are provided.
22//
23// Benchmark
24// ---------
25// PCG: 1000000 Uint64 generated in 15.58s
26// ISAAC: 1000000 Uint64 generated in 8.95s
27// ISAAC: 1000000 Uint32 generated in 7.66s
28// Ratio: x1.74 times faster than PCG (uint64)
29// Ratio: x2.03 times faster than PCG (uint32)
30//
31// Use it directly:
32//
33// prng = isaac.New() // pass 0 to 256 uint64 seeds; if fewer than 256 are provided, the rest
34// // will be generated using the xorshiftr128plus PRNG.
35//
36// Or use it as a drop-in replacement for the default PRNT in Rand:
37//
38// source = isaac64.New()
39// prng := rand.New(source)
40package isaac64
41
42import (
43 "errors"
44 "math"
45
46 "gno.land/p/demo/entropy"
47 "gno.land/p/nt/ufmt"
48 "gno.land/p/wyhaines/rand/xorshiftr128plus"
49)
50
51const (
52 RANDSIZL = 8
53 RANDSIZ = 1 << RANDSIZL // 256
54)
55
56type ISAAC struct {
57 randrsl [256]uint64
58 randcnt uint64
59 mm [256]uint64
60 aa, bb, cc uint64
61 seed [256]uint64
62}
63
64// ISAAC requires a large, 256-element seed. This implementation will leverage the entropy
65// package combined with the xorshiftr128plus PRNG to generate any missing seeds if fewer than
66// the required number of arguments are provided.
67func New(seeds ...uint64) *ISAAC {
68 isaac := &ISAAC{}
69 seed := [256]uint64{}
70
71 index := 0
72 for index = 0; index < len(seeds) && index < 256; index++ {
73 seed[index] = seeds[index]
74 }
75
76 if index < 2 {
77 e := entropy.New()
78 for ; index < 2; index++ {
79 seed[index] = e.Value64()
80 }
81 }
82
83 // Use the first two seeds as seeding inputs for xorshiftr128plus, in order to
84 // use it to provide any remaining missing seeds.
85 prng := xorshiftr128plus.New(
86 seed[0],
87 seed[1],
88 )
89 for ; index < 256; index++ {
90 seed[index] = prng.Uint64()
91 }
92 isaac.Seed(seed)
93 return isaac
94}
95
96// Reinitialize the generator with a new seed. A seed must be composed of 256 uint64.
97func (isaac *ISAAC) Seed(seed [256]uint64) {
98 isaac.randrsl = seed
99 isaac.seed = seed
100 isaac.randinit(true)
101}
102
103// beUint64() decodes a uint64 from a set of eight bytes, assuming big endian encoding.
104func beUint64(b []byte) uint64 {
105 _ = b[7] // bounds check hint to compiler
106 return uint64(b[7]) | uint64(b[6])<<8 | uint64(b[5])<<16 | uint64(b[4])<<24 |
107 uint64(b[3])<<32 | uint64(b[2])<<40 | uint64(b[1])<<48 | uint64(b[0])<<56
108}
109
110// bePutUint64() encodes a uint64 into a buffer of eight bytes.
111func bePutUint64(b []byte, v uint64) {
112 _ = b[7] // early bounds check to guarantee safety of writes below
113 b[0] = byte(v >> 56)
114 b[1] = byte(v >> 48)
115 b[2] = byte(v >> 40)
116 b[3] = byte(v >> 32)
117 b[4] = byte(v >> 24)
118 b[5] = byte(v >> 16)
119 b[6] = byte(v >> 8)
120 b[7] = byte(v)
121}
122
123// A label to identify the marshalled data.
124var marshalISAACLabel = []byte("isaac:")
125
126// MarshalBinary() returns a byte array that encodes the state of the PRNG. This can later be used
127// with UnmarshalBinary() to restore the state of the PRNG.
128// MarshalBinary implements the encoding.BinaryMarshaler interface.
129func (isaac *ISAAC) MarshalBinary() ([]byte, error) {
130 b := make([]byte, 6+2048*3+8*3+8) // 6 + 2048*3 + 8*3 + 8 == 6182
131 copy(b, marshalISAACLabel)
132 offset := 6
133 for i := 0; i < 256; i++ {
134 bePutUint64(b[offset:], isaac.seed[i])
135 offset += 8
136 }
137 for i := 0; i < 256; i++ {
138 bePutUint64(b[offset:], isaac.randrsl[i])
139 offset += 8
140 }
141 for i := 0; i < 256; i++ {
142 bePutUint64(b[offset:], isaac.mm[i])
143 offset += 8
144 }
145 bePutUint64(b[offset:], isaac.aa)
146 offset += 8
147 bePutUint64(b[offset:], isaac.bb)
148 offset += 8
149 bePutUint64(b[offset:], isaac.cc)
150 offset += 8
151 bePutUint64(b[offset:], isaac.randcnt)
152 return b, nil
153}
154
155// errUnmarshalISAAC is returned when unmarshalling fails.
156var errUnmarshalISAAC = errors.New("invalid ISAAC encoding")
157
158// UnmarshalBinary() restores the state of the PRNG from a byte array that was created with MarshalBinary().
159// UnmarshalBinary implements the encoding.BinaryUnmarshaler interface.
160func (isaac *ISAAC) UnmarshalBinary(data []byte) error {
161 if len(data) != 6182 || string(data[:6]) != string(marshalISAACLabel) {
162 return errUnmarshalISAAC
163 }
164 offset := 6
165 for i := 0; i < 256; i++ {
166 isaac.seed[i] = beUint64(data[offset:])
167 offset += 8
168 }
169 for i := 0; i < 256; i++ {
170 isaac.randrsl[i] = beUint64(data[offset:])
171 offset += 8
172 }
173 for i := 0; i < 256; i++ {
174 isaac.mm[i] = beUint64(data[offset:])
175 offset += 8
176 }
177 isaac.aa = beUint64(data[offset:])
178 offset += 8
179 isaac.bb = beUint64(data[offset:])
180 offset += 8
181 isaac.cc = beUint64(data[offset:])
182 offset += 8
183 isaac.randcnt = beUint64(data[offset:])
184 return nil
185}
186
187func (isaac *ISAAC) randinit(flag bool) {
188 var a, b, c, d, e, f, g, h uint64
189 isaac.aa = 0
190 isaac.bb = 0
191 isaac.cc = 0
192
193 a = 0x9e3779b97f4a7c13
194 b = 0x9e3779b97f4a7c13
195 c = 0x9e3779b97f4a7c13
196 d = 0x9e3779b97f4a7c13
197 e = 0x9e3779b97f4a7c13
198 f = 0x9e3779b97f4a7c13
199 g = 0x9e3779b97f4a7c13
200 h = 0x9e3779b97f4a7c13
201
202 // scramble it
203 for i := 0; i < 4; i++ {
204 mix(&a, &b, &c, &d, &e, &f, &g, &h)
205 }
206
207 // fill in mm[] with messy stuff
208 for i := 0; i < RANDSIZ; i += 8 {
209 if flag {
210 a += isaac.randrsl[i]
211 b += isaac.randrsl[i+1]
212 c += isaac.randrsl[i+2]
213 d += isaac.randrsl[i+3]
214 e += isaac.randrsl[i+4]
215 f += isaac.randrsl[i+5]
216 g += isaac.randrsl[i+6]
217 h += isaac.randrsl[i+7]
218 }
219 mix(&a, &b, &c, &d, &e, &f, &g, &h)
220 isaac.mm[i] = a
221 isaac.mm[i+1] = b
222 isaac.mm[i+2] = c
223 isaac.mm[i+3] = d
224 isaac.mm[i+4] = e
225 isaac.mm[i+5] = f
226 isaac.mm[i+6] = g
227 isaac.mm[i+7] = h
228 }
229
230 if flag {
231 // do a second pass to make all of the seed affect all of mm
232 for i := 0; i < RANDSIZ; i += 8 {
233 a += isaac.mm[i]
234 b += isaac.mm[i+1]
235 c += isaac.mm[i+2]
236 d += isaac.mm[i+3]
237 e += isaac.mm[i+4]
238 f += isaac.mm[i+5]
239 g += isaac.mm[i+6]
240 h += isaac.mm[i+7]
241 mix(&a, &b, &c, &d, &e, &f, &g, &h)
242 isaac.mm[i] = a
243 isaac.mm[i+1] = b
244 isaac.mm[i+2] = c
245 isaac.mm[i+3] = d
246 isaac.mm[i+4] = e
247 isaac.mm[i+5] = f
248 isaac.mm[i+6] = g
249 isaac.mm[i+7] = h
250 }
251 }
252
253 isaac.isaac()
254 isaac.randcnt = RANDSIZ
255}
256
257func mix(a, b, c, d, e, f, g, h *uint64) {
258 *a -= *e
259 *f ^= *h >> 9
260 *h += *a
261
262 *b -= *f
263 *g ^= *a << 9
264 *a += *b
265
266 *c -= *g
267 *h ^= *b >> 23
268 *b += *c
269
270 *d -= *h
271 *a ^= *c << 15
272 *c += *d
273
274 *e -= *a
275 *b ^= *d >> 14
276 *d += *e
277
278 *f -= *b
279 *c ^= *e << 20
280 *e += *f
281
282 *g -= *c
283 *d ^= *f >> 17
284 *f += *g
285
286 *h -= *d
287 *e ^= *g << 14
288 *g += *h
289}
290
291func ind(mm []uint64, x uint64) uint64 {
292 return mm[(x>>3)&(RANDSIZ-1)]
293}
294
295func (isaac *ISAAC) isaac() {
296 var a, b, x, y uint64
297 a = isaac.aa
298 b = isaac.bb + isaac.cc + 1
299 isaac.cc++
300
301 m := isaac.mm[:]
302 r := isaac.randrsl[:]
303
304 var i, m2Index int
305
306 // First half
307 for i = 0; i < RANDSIZ/2; i++ {
308 m2Index = i + RANDSIZ/2
309 switch i % 4 {
310 case 0:
311 a = ^(a ^ (a << 21)) + m[m2Index]
312 case 1:
313 a = (a ^ (a >> 5)) + m[m2Index]
314 case 2:
315 a = (a ^ (a << 12)) + m[m2Index]
316 case 3:
317 a = (a ^ (a >> 33)) + m[m2Index]
318 }
319 x = m[i]
320 y = ind(m, x) + a + b
321 m[i] = y
322 b = ind(m, y>>RANDSIZL) + x
323 r[i] = b
324 }
325
326 // Second half
327 for i = RANDSIZ / 2; i < RANDSIZ; i++ {
328 m2Index = i - RANDSIZ/2
329 switch i % 4 {
330 case 0:
331 a = ^(a ^ (a << 21)) + m[m2Index]
332 case 1:
333 a = (a ^ (a >> 5)) + m[m2Index]
334 case 2:
335 a = (a ^ (a << 12)) + m[m2Index]
336 case 3:
337 a = (a ^ (a >> 33)) + m[m2Index]
338 }
339 x = m[i]
340 y = ind(m, x) + a + b
341 m[i] = y
342 b = ind(m, y>>RANDSIZL) + x
343 r[i] = b
344 }
345
346 isaac.bb = b
347 isaac.aa = a
348}
349
350// Return a 64 bit random integer.
351func (isaac *ISAAC) Uint64() uint64 {
352 if isaac.randcnt == 0 {
353 isaac.isaac()
354 isaac.randcnt = RANDSIZ
355 }
356 isaac.randcnt--
357 return isaac.randrsl[isaac.randcnt]
358}
359
360var gencycle int = 0
361var bufferFor32 uint64 = uint64(0)
362
363// Return a 32 bit random integer, composed of the high 32 bits of the generated 32 bit result.
364func (isaac *ISAAC) Uint32() uint32 {
365 if gencycle == 0 {
366 bufferFor32 = isaac.Uint64()
367 gencycle = 1
368 return uint32(bufferFor32 >> 32)
369 }
370
371 gencycle = 0
372 return uint32(bufferFor32 & 0xffffffff)
373}
374
375// Until there is better benchmarking support in gno, you can test the performance of this PRNG with this function.
376// This isn't perfect, since it will include the startup time of gno in the results, but this will give you a timing
377// for generating a million random uint64 numbers on any unix based system:
378//
379// `time gno run -expr 'benchmarkISAAC()' isaac64.gno
380func benchmarkISAAC(_iterations ...int) {
381 iterations := 1000000
382 if len(_iterations) > 0 {
383 iterations = _iterations[0]
384 }
385 isaac := New()
386
387 for i := 0; i < iterations; i++ {
388 _ = isaac.Uint64()
389 }
390 ufmt.Println(ufmt.Sprintf("ISAAC: generated %d uint64\n", iterations))
391}
392
393// The averageISAAC() function is a simple benchmarking helper to demonstrate
394// the most basic statistical property of the ISAAC PRNG.
395func averageISAAC(_iterations ...int) {
396 target := uint64(500000)
397 iterations := 1000000
398
399 ufmt.Println(
400 ufmt.Sprintf(
401 "Averaging %d random numbers. The average should be very close to %d.\n",
402 iterations,
403 target))
404
405 if len(_iterations) > 0 {
406 iterations = _iterations[0]
407 }
408 isaac := New(987654321987654321, 123456789987654321, 1, 997755331886644220)
409
410 var average float64 = 0
411 var squares []uint64 = make([]uint64, iterations)
412 for i := 0; i < iterations; i++ {
413 n := isaac.Uint64()%(target*2) + 1
414 average += (float64(n) - average) / float64(i+1)
415 squares[i] = n
416 }
417
418 sum_of_squares := uint64(0)
419 // transform numbers into their squares of the distance from the average
420 for i := 0; i < iterations; i++ {
421 difference := average - float64(squares[i])
422 square := uint64(difference * difference)
423 sum_of_squares += square
424 }
425
426 ufmt.Println(ufmt.Sprintf("ISAAC average of %d uint64: %f\n", iterations, average))
427 ufmt.Println(ufmt.Sprintf("ISAAC standard deviation : %f\n", math.Sqrt(float64(sum_of_squares)/float64(iterations))))
428 ufmt.Println(ufmt.Sprintf("ISAAC theoretical perfect deviation: %f\n", (float64(target*2)-1)/math.Sqrt(12)))
429}
