1	import Base
2	import ../../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	
22	def 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
30	def 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	
37	def lcg(+mul: Nat, +inc: Nat, +s: Nat) -> Nat:
38	  C.low(128n, Nat.add(Nat.mul(s, mul), inc))
39	
40	def 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
44	def dxsm1(+cm: Nat, +hi: Nat) -> Nat:
45	  C.low(64n, Nat.mul(xor_bits(64n, hi, C.high(32n, hi)), cm))
46	
47	def dxsm(+cm: Nat, +hi: Nat, +lo: Nat) -> Nat:
48	  dxsm3(dxsm1(cm, hi), lo)