1	import Base
2	import ../../../spec/lib/common.bend as C
3	import ../../../spec/math/w64.bend as SW
4	import ../../../spec/math/random/pcg.bend as SP
5	import ../../../spec/math/random.bend as SRM
6	import ../../../src/math/u64.bend as WU
7	import ../../../src/math/w64.bend as X
8	import ../../../src/math/random/pcg.bend as P
9	import ../../lib/lemmas/proofs/nat_algebra.bend as NA
10	import ../typed/width.bend as WW
11	import ../typed/w64add.bend as WA
12	import ../typed/w64sh.bend as SH
13	import ./pcg/xor.bend as XO
14	import ../../lib/lemmas/spec/numeric.bend as S
15	import ../typed/w64m128.bend as M128
16	import ./pcg/step.bend as ST
17	import ./pcg/impl.bend as IM
18	
19	# Entry point: `bend proofs/math/random/proof_pcg.bend` checks the PCG clauses
20	# of spec/math/random.bend under the clause's name, for every state and every
21	# multiplier and increment; no holes, no axioms. One step of
22	# src/math/random/pcg.bend is the 128-bit LCG (the width-generic algebra is in
23	# pcg/step.bend, instantiated at 64 in pcg/impl.bend) and its DXSM output is
24	# spec/math/random/pcg.bend's dxsm: xor with a right shift, products mod 2^64
25	# and "or 1" on the low word, each by its value lemma. The factors of each
26	# product stay variables (mulg), so no product is ever normalized, and the
27	# clause types are stated once (a clause type is expensive to normalize).
28	# The other clauses are checked by proofs/math/random/proof.bend,
29	# proof_draws.bend and proof_float.bend.
30	
31	def PCG.step(+mh: WU.U64, +ml: WU.U64, +ih: WU.U64, +il: WU.U64, +p: P.PCG) -> SRM.PCG.step(mh, ml, ih, il, p):
32	  match p:
33	    case P.P{+hi, +lo}:
34	      IM.step_pair(mh, ml, ih, il, hi, lo, X.mul128(lo, ml), IM.m128(lo, ml))
35	
36	def PCG.constants(+p: P.PCG) -> SRM.PCG.constants(p):
37	  {==}
38	
39	def v(+x: U32) -> Nat:
40	  U32.to_nat(x)
41	
42	# x ^ (x >> k)
43	def xs(+a: WU.U64, +k: Nat) -> {SW.value(P.xor64(a, X.shr(a, k))) == SP.xor_bits(64n, SW.value(a), C.high(k, SW.value(a))) : Nat}:
44	  Equal.trans(Nat, SW.value(P.xor64(a, X.shr(a, k))), SP.xor_bits(64n, SW.value(a), SW.value(X.shr(a, k))), SP.xor_bits(64n, SW.value(a), C.high(k, SW.value(a))),
45	    XO.xor64_value(a, X.shr(a, k)),
46	    Equal.cong(Nat, Nat, z => SP.xor_bits(64n, SW.value(a), z), SW.value(X.shr(a, k)), C.high(k, SW.value(a)), SH.shr_value(a, k)))
47	
48	def plus1(+x: Nat) -> {Nat.add(x, 1n) == 1n+x : Nat}:
49	  Equal.trans(Nat, Nat.add(x, 1n), 1n+Nat.add(x, 0n), 1n+x, NA.add_succ(x, 0n), Equal.cong(Nat, Nat, z => 1n+z, Nat.add(x, 0n), x, NA.add_zero(x)))
50	
51	# lo | 1
52	def or_value(+lo: WU.U64) -> {SW.value(WU.U64{U32.or(X.lo(lo), 1), X.hi(lo)}) == Nat.add(Nat.double(C.half(SW.value(lo))), 1n) : Nat}:
53	  match lo:
54	    case WU.U64{+l, +h}:
55	      +vl = v(l)
56	      +vh = v(h)
57	      +D = Nat.add(Nat.double(C.half(vl)), C.shift(32n, vh))
58	      +e1 = Equal.cong(Nat, Nat, z => Nat.add(z, C.shift(32n, vh)), v(U32.or(l, 1)), 1n+Nat.double(C.half(vl)), SH.or1(l))
59	      +e2 = Equal.cong(Nat, Nat, z => 1n+z, D, Nat.double(Nat.add(C.half(vl), C.shift(31n, vh))), Equal.sym(Nat, Nat.double(Nat.add(C.half(vl), C.shift(31n, vh))), D, NA.double_add(C.half(vl), C.shift(31n, vh))))
60	      +e3 = Equal.cong(Nat, Nat, z => 1n+Nat.double(z), Nat.add(C.half(vl), C.shift(31n, vh)), C.half(Nat.add(vl, C.shift(32n, vh))), Equal.sym(Nat, C.half(Nat.add(vl, C.shift(32n, vh))), Nat.add(C.half(vl), C.shift(31n, vh)), WW.half_dbl(vl, C.shift(31n, vh))))
61	      +e4 = Equal.sym(Nat, Nat.add(Nat.double(C.half(Nat.add(vl, C.shift(32n, vh)))), 1n), 1n+Nat.double(C.half(Nat.add(vl, C.shift(32n, vh)))), plus1(Nat.double(C.half(Nat.add(vl, C.shift(32n, vh))))))
62	      Equal.trans(Nat, Nat.add(v(U32.or(l, 1)), C.shift(32n, vh)), 1n+D, Nat.add(Nat.double(C.half(Nat.add(vl, C.shift(32n, vh)))), 1n), e1,
63	        Equal.trans(Nat, 1n+D, 1n+Nat.double(Nat.add(C.half(vl), C.shift(31n, vh))), Nat.add(Nat.double(C.half(Nat.add(vl, C.shift(32n, vh)))), 1n), e2,
64	          Equal.trans(Nat, 1n+Nat.double(Nat.add(C.half(vl), C.shift(31n, vh))), 1n+Nat.double(C.half(Nat.add(vl, C.shift(32n, vh)))), Nat.add(Nat.double(C.half(Nat.add(vl, C.shift(32n, vh)))), 1n), e3, e4)))
65	
66	# the value of a product mod 2^64 from the values of its factors; the
67	# factors stay variables here, so no product is ever normalized
68	def mulg(+x: WU.U64, +o: WU.U64, +xv: Nat, +hx: {SW.value(x) == xv : Nat}, +ov: Nat, +ho: {SW.value(o) == ov : Nat}) -> {SW.value(X.mul(x, o)) == C.low(64n, Nat.mul(xv, ov)) : Nat}:
69	  Equal.trans(Nat, SW.value(X.mul(x, o)), C.low(64n, Nat.mul(SW.value(x), SW.value(o))), C.low(64n, Nat.mul(xv, ov)), WA.mul_value(x, o),
70	    Equal.trans(Nat, C.low(64n, Nat.mul(SW.value(x), SW.value(o))), C.low(64n, Nat.mul(xv, SW.value(o))), C.low(64n, Nat.mul(xv, ov)),
71	      Equal.cong(Nat, Nat, z => C.low(64n, Nat.mul(z, SW.value(o))), SW.value(x), xv, hx),
72	      Equal.cong(Nat, Nat, z => C.low(64n, Nat.mul(xv, z)), SW.value(o), ov, ho)))
73	
74	def xs32(+h: WU.U64) -> {SW.value(P.xor64(h, X.shr(h, 32n))) == SP.xor_bits(64n, SW.value(h), C.high(32n, SW.value(h))) : Nat}:
75	  xs(h, 32n)
76	
77	def xs48(+h: WU.U64) -> {SW.value(P.xor64(h, X.shr(h, 48n))) == SP.xor_bits(64n, SW.value(h), C.high(48n, SW.value(h))) : Nat}:
78	  xs(h, 48n)
79	
80	# x ^ (x >> 48) for a word of value t
81	def xt(+H: WU.U64, +t: Nat, +ht: {SW.value(H) == t : Nat}) -> {SW.value(P.xor64(H, X.shr(H, 48n))) == SP.xor_bits(64n, t, C.high(48n, t)) : Nat}:
82	  Equal.trans(Nat, SW.value(P.xor64(H, X.shr(H, 48n))), SP.xor_bits(64n, SW.value(H), C.high(48n, SW.value(H))), SP.xor_bits(64n, t, C.high(48n, t)),
83	    xs48(H), Equal.cong(Nat, Nat, z => SP.xor_bits(64n, z, C.high(48n, z)), SW.value(H), t, ht))
84	
85	# the first half's value is t
86	def half1(+cm: WU.U64, +hi: WU.U64, +t: Nat, +ht: {t == SP.dxsm1(SW.value(cm), SW.value(hi)) : Nat}) -> {SW.value(X.mul(P.xor64(hi, X.shr(hi, 32n)), cm)) == t : Nat}:
87	  Equal.trans(Nat, SW.value(X.mul(P.xor64(hi, X.shr(hi, 32n)), cm)), SP.dxsm1(SW.value(cm), SW.value(hi)), t,
88	    mulg(P.xor64(hi, X.shr(hi, 32n)), cm, SP.xor_bits(64n, SW.value(hi), C.high(32n, SW.value(hi))), xs32(hi), SW.value(cm), {==}),
89	    Equal.sym(Nat, t, SP.dxsm1(SW.value(cm), SW.value(hi)), ht))
90	
91	# THEOREM (PCG.output)
92	# values are known (x ^ x >> 48 with x of value t, and lo | 1)
93	def PCG.output(+cm: WU.U64, +hi: WU.U64, +lo: WU.U64, +t: Nat, +ht: {t == SP.dxsm1(SW.value(cm), SW.value(hi)) : Nat}) -> SRM.PCG.output(cm, hi, lo, t, ht):
94	  +H = X.mul(P.xor64(hi, X.shr(hi, 32n)), cm)
95	  mulg(P.xor64(H, X.shr(H, 48n)), WU.U64{U32.or(X.lo(lo), 1), X.hi(lo)}, SP.xor_bits(64n, t, C.high(48n, t)), xt(H, t, half1(cm, hi, t, ht)), Nat.add(Nat.double(C.half(SW.value(lo))), 1n), or_value(lo))