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