1-- A deterministic generator of pseudorandom numbers: xoshiro256** 1.0,
2-- by David Blackman and Sebastiano Vigna, "Scrambled linear pseudorandom
3-- number generators", ACM Transactions on Mathematical Software 47(4),
4-- 2021, from their public-domain reference implementation at
5-- https://prng.di.unimi.it/xoshiro256starstar.c, seeded from a 64-bit
6-- seed through their splitmix64.c, as they suggest. It is written with
7-- shifts, exclusive ors and additions on `u64`, so a seed gives the same
8-- sequence on every target, a 32-bit Cortex-M0 included, and its
9-- multiplications by 5 and 9 are a shift and an addition each.
10
11import core/io
12
13state_value: type = struct
14 s0: u64
15 s1: u64
16 s2: u64
17 s3: u64
18end state_value
19
20--- A generator: 256 bits of state that each draw advances, a value the
21--- caller keeps and passes `inout`. A seed always gives the same
22--- sequence, on every target, which is what makes it useful for tests,
23--- simulations and sampling; a generator is not cryptographic, and its
24--- outputs reveal its state, so never make a key, a token or a password
25--- with one: take those from `io.fill_random` directly. Copying a
26--- generator copies its state, and the copy repeats what the original
27--- will draw.
28public generator: type = state_value
29
30-- splitmix64's step: advance `seed` by the golden ratio and mix it.
31splitmix: (inout seed: u64) -> (value: u64) =
32 seed = seed +% 0x9e3779b97f4a7c15
33 mut z: u64 = seed
34 z = (z ^ (z >> 30)) *% 0xbf58476d1ce4e5b9
35 z = (z ^ (z >> 27)) *% 0x94d049bb133111eb
36 value = z ^ (z >> 31)
37end splitmix
38
39--- A generator seeded from `seed`: its state is four outputs of
40--- splitmix64 started at `seed`, as the authors suggest, which are never
41--- all zero. Seed 0 gives first 11091344671253066420, as the reference
42--- implementation does.
43public new: (seed: u64) -> (result: generator) =
44 mut at: u64 = seed
45 first: u64 = splitmix(at)
46 second: u64 = splitmix(at)
47 third: u64 = splitmix(at)
48 fourth: u64 = splitmix(at)
49 result = (s0: first, s1: second, s2: third, s3: fourth)
50end new
51
52--- A generator seeded from 32 bytes of the provider's entropy, read as
53--- four little-endian `u64`s, so no two runs repeat each other. The one
54--- state xoshiro256** cannot leave, all zero, is replaced by `new(0)`'s.
55--- It fails as the provider fails.
56public from_entropy: (provider: type is io.entropy, host: provider)
57 -> (result: generator) ! io.io_failed =
58 mut bytes: [32]u8 = zeroed
59 try io.fill_random(host, bytes[0..<32])
60 result = (s0: little_endian(bytes[0..<8]),
61 s1: little_endian(bytes[8..<16]),
62 s2: little_endian(bytes[16..<24]),
63 s3: little_endian(bytes[24..<32]))
64 if result.s0 == 0 and result.s1 == 0 and result.s2 == 0
65 and result.s3 == 0 then
66 result = new(0)
67 end if
68end from_entropy
69
70-- The eight bytes of `bytes`, least significant first.
71little_endian: (bytes: []u8) -> (value: u64) =
72 value = 0
73 for index in usize(0)..<8 do
74 value = value | (u64(bytes[index]) << u64(8 * index))
75 end for
76end little_endian
77
78-- `value` rotated left by `count` bits, 0 < count < 64.
79rotate: (value: u64, count: u64) -> (result: u64) =
80 result = (value << count) | (value >> (64 - count))
81end rotate
82
83--- The next 64 bits of the sequence, each value equally likely.
84public next: (inout source: generator) -> (value: u64) =
85 s1: u64 = source.s1
86 times_five: u64 = (s1 << 2) +% s1
87 rotated: u64 = rotate(times_five, 7)
88 value = (rotated << 3) +% rotated
89 t: u64 = s1 << 17
90 source.s2 = source.s2 ^ source.s0
91 source.s3 = source.s3 ^ source.s1
92 source.s1 = source.s1 ^ source.s2
93 source.s0 = source.s0 ^ source.s3
94 source.s2 = source.s2 ^ t
95 source.s3 = rotate(source.s3, 45)
96end next
97
98--- The next 32 bits of the sequence: the upper half of `next`, whose bits
99--- the authors recommend over the lower.
100public next_u32: (inout source: generator) -> (value: u32) =
101 value = u32(next(source) >> 32)
102end next_u32
103
104--- A value from 0 to `bound - 1`, each equally likely: the next draw's
105--- bits below `bound`'s highest bit, drawn again while they are `bound`
106--- or more, which is never more often than not, so no value is favoured
107--- as a remainder would favour some. It uses no division. A `bound` of
108--- zero has no value below it and traps.
109public below: (inout source: generator, bound: u64) -> (value: u64) =
110 if bound == 0 then
111 compiler.unreachable()
112 end if
113 mut mask: u64 = bound - 1
114 mask = mask | (mask >> 1)
115 mask = mask | (mask >> 2)
116 mask = mask | (mask >> 4)
117 mask = mask | (mask >> 8)
118 mask = mask | (mask >> 16)
119 mask = mask | (mask >> 32)
120 loop do
121 value = next(source) & mask
122 return when value < bound
123 end loop
124end below
125
126--- Put `values` in an order drawn uniformly from every order they could
127--- have, by Fisher and Yates's shuffle: from the last, each is swapped
128--- with one drawn by `below` from those before it and itself.
129public shuffle: (item: type, inout source: generator, values: []mut item)
130 -> none =
131 mut at: usize = lenof values
132 while at > 1 do
133 chosen: usize = usize(below(source, u64(at)))
134 dec at
135 saved: item = values[at]
136 values[at] = values[chosen]
137 values[chosen] = saved
138 end while
139end shuffle