spec/math/random/pcg.bend
1import Base
2import bend-collections-laws-math@1.0.0.0/spec/lib/common.bend as C
3
4# Executable specification of Go's PCG (math/rand/v2 pcg.go; O'Neill, "PCG:
5# A Family of Simple Fast Space-Efficient Statistically Good Algorithms for
6# Random Number Generation", 2014, with the DXSM output function of numpy
7# and Go) on natural numbers:
8#
9# lcg(mul, inc, s) one step of the 128-bit state: s * mul + inc
10# mod 2^128
11# dxsm(cm, hi, lo) the output of the state hi * 2^64 + lo:
12# hi ^= hi >> 32; hi *= cm (dxsm1); then
13# hi ^= hi >> 48; hi *= lo | 1 (dxsm3), each
14# product mod 2^64
15#
16# Go's constants are mul = 2549297995355413924 * 2^64 + 4865540595714422341,
17# inc = 6364136223846793005 * 2^64 + 1442695040888963407 and
18# cm = 0xda942042e4dd58b5 (src/math/random/pcg.bend holds them as 64-bit
19# words); the functions take them as parameters, so no closed 64-bit
20# constant is ever expanded by the checker (spec/lib/common.bend).
21
22def b2n(b: Bool) -> Nat:
23 match b:
24 case True{}:
25 1n
26 case False{}:
27 0n
28
29# the exclusive or of the low w bits of a and b
30def xor_bits(w: Nat, +a: Nat, +b: Nat) -> Nat:
31 match w:
32 case 0n:
33 0n
34 case 1n+k:
35 Nat.add(b2n(Bool.not(Nat.is_eq(C.bit(a), C.bit(b)))), Nat.double(xor_bits(k, C.half(a), C.half(b))))
36
37def lcg(+mul: Nat, +inc: Nat, +s: Nat) -> Nat:
38 C.low(128n, Nat.add(Nat.mul(s, mul), inc))
39
40def dxsm3(+h: Nat, +lo: Nat) -> Nat:
41 C.low(64n, Nat.mul(xor_bits(64n, h, C.high(48n, h)), Nat.add(Nat.double(C.half(lo)), 1n)))
42
43# the first half: hi ^= hi >> 32; hi *= cm
44def dxsm1(+cm: Nat, +hi: Nat) -> Nat:
45 C.low(64n, Nat.mul(xor_bits(64n, hi, C.high(32n, hi)), cm))
46
47def dxsm(+cm: Nat, +hi: Nat, +lo: Nat) -> Nat:
48 dxsm3(dxsm1(cm, hi), lo)