proofs/math/random/proof_pcg.bend

bend-collections-laws-math@1.0.0.0 · 95 lines · 12 definitions

1import Base
2import bend-collections-laws-math@1.0.0.0/spec/lib/common.bend as C
3import bend-collections-laws-math@1.0.0.0/spec/math/w64.bend as SW
4import bend-collections-laws-math@1.0.0.0/spec/math/random/pcg.bend as SP
5import bend-collections-laws-math@1.0.0.0/spec/math/random.bend as SRM
6import bend-collections-laws-math@1.0.0.0/src/math/u64.bend as WU
7import bend-collections-laws-math@1.0.0.0/src/math/w64.bend as X
8import bend-collections-laws-math@1.0.0.0/src/math/random/pcg.bend as P
9import bend-collections-laws-math@1.0.0.0/proofs/lib/lemmas/proofs/nat_algebra.bend as NA
10import bend-collections-laws-math@1.0.0.0/proofs/math/typed/width.bend as WW
11import bend-collections-laws-math@1.0.0.0/proofs/math/typed/w64add.bend as WA
12import bend-collections-laws-math@1.0.0.0/proofs/math/typed/w64sh.bend as SH
13import bend-collections-laws-math@1.0.0.0/proofs/math/random/pcg/xor.bend as XO
14import bend-collections-laws-math@1.0.0.0/proofs/lib/lemmas/spec/numeric.bend as S
15import bend-collections-laws-math@1.0.0.0/proofs/math/typed/w64m128.bend as M128
16import bend-collections-laws-math@1.0.0.0/proofs/math/random/pcg/step.bend as ST
17import bend-collections-laws-math@1.0.0.0/proofs/math/random/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
31def 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
36def PCG.constants(+p: P.PCG) -> SRM.PCG.constants(p):
37  {==}
38
39def v(+x: U32) -> Nat:
40  U32.to_nat(x)
41
42# x ^ (x >> k)
43def 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
48def 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
52def 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
68def 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
74def 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
77def 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
81def 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
86def 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)
93def 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))