1	import Base
2	import ../../../spec/lib/common.bend as C
3	import ../../../spec/math/random/rand.bend as SR
4	import ../../lib/logic.bend as L
5	import ../../lib/nat.bend as N
6	import ../typed/width.bend as WI
7	import ../natural/arith.bend as AR
8	import ../../lib/arith.bend as LA
9	import ../../lib/u32alg.bend as UA
10	import ../../lib/lemmas/proofs/nat_algebra.bend as NA
11	
12	# Lemire's rejection is exactly unbiased (Lemire, "Fast Random Integer
13	# Generation in an Interval", ACM TOMACS 2019, section 4): for every width w,
14	# bound 0 < n < 2^w and k < n, exactly floor(2^w / n) of the 2^w source
15	# outputs x draw k under Go's uint64n decision (spec/math/random/rand.bend
16	# draw), in both of its branches.
17	#
18	# Lemire's branch: x draws k exactly when k 2^w + t <= x n < (k + 1) 2^w,
19	# t = 2^w mod n. Counting x < N with x n < c gives min(N, ceil(c / n)), and
20	# (k + 1) 2^w - (k 2^w + t) = 2^w - t = n floor(2^w / n), so the two ceilings
21	# differ by floor(2^w / n).
22	# The mask branch (n = 2^j, Go's n & (n - 1) == 0): x & (n - 1) is the low j
23	# bits of x, and splitting x = 2y + b halves the width and the bound.
24	
25	def b2n(b: Bool) -> Nat:
26	  SR.b2n(b)
27	
28	# ---- Bool helpers ----
29	
30	def and_false(+b: Bool) -> {Bool.and(b, False{}) == False{} : Bool}:
31	  match b:
32	    case True{}:
33	      {==}
34	    case False{}:
35	      {==}
36	
37	# lo < n and lo < t is lo < t, when t < n
38	def and_lt(+lo: Nat, +n: Nat, +t: Nat, +ht: {Nat.is_lt(t, n) == True{} : Bool}, +c: Bool, +hc: {Nat.is_lt(lo, t) == c : Bool}) -> {Bool.and(Nat.is_lt(lo, n), Nat.is_lt(lo, t)) == c : Bool}:
39	  match c:
40	    case True{}:
41	      %Equal.sym(Bool, Nat.is_lt(lo, n), True{}, N.lt_trans(lo, t, n, hc, ht)) : {Bool.and(_, Nat.is_lt(lo, t)) == True{} : Bool}
42	      hc
43	    case False{}:
44	      %Equal.sym(Bool, Nat.is_lt(lo, t), False{}, hc) : {Bool.and(Nat.is_lt(lo, n), _) == False{} : Bool}
45	      and_false(Nat.is_lt(lo, n))
46	
47	# the hit of an accepted or rejected hi
48	def hit_acc(+hi: Nat, +k: Nat, +ok: Bool) -> Nat:
49	  SR.hit(SR.accept(hi, ok), k)
50	
51	# s <= t + s
52	def le_plus(+t: Nat, +s: Nat) -> {Nat.is_le(s, Nat.add(t, s)) == True{} : Bool}:
53	  %N.add_comm(s, t) : {Nat.is_le(s, _) == True{} : Bool}
54	  N.le_add_right(s, t)
55	
56	# ---- one output ----
57	
58	def fin_lt(+hi: Nat, +k: Nat, +ok: Bool, +h: {Nat.is_lt(hi, k) == True{} : Bool}) -> {Nat.add(hit_acc(hi, k, ok), 1n) == 1n : Nat}:
59	  match ok:
60	    case True{}:
61	      %Equal.sym(Bool, Nat.is_eq(hi, k), False{}, N.is_eq_lt(hi, k, h)) : {Nat.add(b2n(_), 1n) == 1n : Nat}
62	      {==}
63	    case False{}:
64	      {==}
65	
66	# hi < k: rejected or not, x n lies below both thresholds
67	def term_lt(+w: Nat, +k: Nat, +t: Nat, +lo: Nat, +hi: Nat, +ok: Bool, +hlo: {C.fits(w, lo) == True{} : Bool}, +hab: {Nat.is_le(Nat.add(t, C.shift(w, k)), C.shift(w, 1n+k)) == True{} : Bool}, +h: {Nat.is_lt(hi, k) == True{} : Bool}) -> {Nat.add(hit_acc(hi, k, ok), b2n(Nat.is_lt(Nat.add(lo, C.shift(w, hi)), Nat.add(t, C.shift(w, k))))) == b2n(Nat.is_lt(Nat.add(lo, C.shift(w, hi)), C.shift(w, 1n+k))) : Nat}:
68	  +y = Nat.add(lo, C.shift(w, hi))
69	  +a = Nat.add(t, C.shift(w, k))
70	  +ya = N.lt_le_trans(y, C.shift(w, 1n+hi), a, WI.lt_hi(w, lo, hi, 0n, 1n+hi, hlo, N.lt_succ(hi)), N.le_trans(C.shift(w, 1n+hi), C.shift(w, k), a, WI.shift_mono(w, 1n+hi, k, N.lt_succ_le_succ(hi, k, h)), le_plus(t, C.shift(w, k))))
71	  +yb = N.lt_le_trans(y, a, C.shift(w, 1n+k), ya, hab)
72	  %Equal.sym(Bool, Nat.is_lt(y, a), True{}, ya) : {Nat.add(hit_acc(hi, k, ok), b2n(_)) == b2n(Nat.is_lt(y, C.shift(w, 1n+k))) : Nat}
73	  %Equal.sym(Bool, Nat.is_lt(y, C.shift(w, 1n+k)), True{}, yb) : {Nat.add(hit_acc(hi, k, ok), 1n) == b2n(_) : Nat}
74	  fin_lt(hi, k, ok, h)
75	
76	def fin_eq(+k: Nat, +c: Bool) -> {Nat.add(hit_acc(k, k, Bool.not(c)), b2n(c)) == 1n : Nat}:
77	  match c:
78	    case True{}:
79	      {==}
80	    case False{}:
81	      %Equal.sym(Bool, Nat.is_eq(k, k), True{}, N.is_eq_refl(k)) : {Nat.add(b2n(_), 0n) == 1n : Nat}
82	      {==}
83	
84	# hi == k: x n is below the upper threshold, and below the lower one exactly
85	# when it is rejected
86	def term_eq(+w: Nat, +k: Nat, +t: Nat, +lo: Nat, +hlo: {C.fits(w, lo) == True{} : Bool}) -> {Nat.add(hit_acc(k, k, Bool.not(Nat.is_lt(lo, t))), b2n(Nat.is_lt(Nat.add(lo, C.shift(w, k)), Nat.add(t, C.shift(w, k))))) == b2n(Nat.is_lt(Nat.add(lo, C.shift(w, k)), C.shift(w, 1n+k))) : Nat}:
87	  %Equal.sym(Bool, Nat.is_lt(Nat.add(lo, C.shift(w, k)), Nat.add(t, C.shift(w, k))), Nat.is_lt(lo, t), WI.lt_cancel_r(lo, t, C.shift(w, k))) : {Nat.add(hit_acc(k, k, Bool.not(Nat.is_lt(lo, t))), b2n(_)) == b2n(Nat.is_lt(Nat.add(lo, C.shift(w, k)), C.shift(w, 1n+k))) : Nat}
88	  %Equal.sym(Bool, Nat.is_lt(Nat.add(lo, C.shift(w, k)), C.shift(w, 1n+k)), True{}, WI.lt_hi(w, lo, k, 0n, 1n+k, hlo, N.lt_succ(k))) : {Nat.add(hit_acc(k, k, Bool.not(Nat.is_lt(lo, t))), b2n(Nat.is_lt(lo, t))) == b2n(_) : Nat}
89	  fin_eq(k, Nat.is_lt(lo, t))
90	
91	def fin_gt(+hi: Nat, +k: Nat, +ok: Bool, +h: {Nat.is_lt(k, hi) == True{} : Bool}) -> {Nat.add(hit_acc(hi, k, ok), 0n) == 0n : Nat}:
92	  match ok:
93	    case True{}:
94	      %Equal.sym(Bool, Nat.is_eq(hi, k), False{}, N.is_eq_sym_false(k, hi, N.is_eq_lt(k, hi, h))) : {Nat.add(b2n(_), 0n) == 0n : Nat}
95	      {==}
96	    case False{}:
97	      {==}
98	
99	# hi > k: x n is at or above both thresholds
100	def term_gt(+w: Nat, +k: Nat, +t: Nat, +lo: Nat, +hi: Nat, +ok: Bool, +hab: {Nat.is_le(Nat.add(t, C.shift(w, k)), C.shift(w, 1n+k)) == True{} : Bool}, +h: {Nat.is_lt(k, hi) == True{} : Bool}) -> {Nat.add(hit_acc(hi, k, ok), b2n(Nat.is_lt(Nat.add(lo, C.shift(w, hi)), Nat.add(t, C.shift(w, k))))) == b2n(Nat.is_lt(Nat.add(lo, C.shift(w, hi)), C.shift(w, 1n+k))) : Nat}:
101	  +y = Nat.add(lo, C.shift(w, hi))
102	  +a = Nat.add(t, C.shift(w, k))
103	  +by = N.le_trans(C.shift(w, 1n+k), C.shift(w, hi), y, WI.shift_mono(w, 1n+k, hi, N.lt_succ_le_succ(k, hi, h)), le_plus(lo, C.shift(w, hi)))
104	  %Equal.sym(Bool, Nat.is_lt(y, a), False{}, N.le_not_lt(y, a, N.le_trans(a, C.shift(w, 1n+k), y, hab, by))) : {Nat.add(hit_acc(hi, k, ok), b2n(_)) == b2n(Nat.is_lt(y, C.shift(w, 1n+k))) : Nat}
105	  %Equal.sym(Bool, Nat.is_lt(y, C.shift(w, 1n+k)), False{}, N.le_not_lt(y, C.shift(w, 1n+k), by)) : {Nat.add(hit_acc(hi, k, ok), 0n) == b2n(_) : Nat}
106	  fin_gt(hi, k, ok, h)
107	
108	def term_cases(+w: Nat, +n: Nat, +k: Nat, +t: Nat, +lo: Nat, +hi: Nat, +hlo: {C.fits(w, lo) == True{} : Bool}, +ht: {Nat.is_lt(t, n) == True{} : Bool}, +hab: {Nat.is_le(Nat.add(t, C.shift(w, k)), C.shift(w, 1n+k)) == True{} : Bool}, +lt: Bool, +eq: Bool, +hlt: {Nat.is_lt(hi, k) == lt : Bool}, +heq: {Nat.is_eq(hi, k) == eq : Bool}) -> {Nat.add(hit_acc(hi, k, Bool.not(Bool.and(Nat.is_lt(lo, n), Nat.is_lt(lo, t)))), b2n(Nat.is_lt(Nat.add(lo, C.shift(w, hi)), Nat.add(t, C.shift(w, k))))) == b2n(Nat.is_lt(Nat.add(lo, C.shift(w, hi)), C.shift(w, 1n+k))) : Nat}:
109	  match lt eq:
110	    case True{} _:
111	      term_lt(w, k, t, lo, hi, Bool.not(Bool.and(Nat.is_lt(lo, n), Nat.is_lt(lo, t))), hlo, hab, hlt)
112	    case False{} True{}:
113	      %Equal.sym(Nat, hi, k, N.eq_from_is_eq(hi, k, heq)) : {Nat.add(hit_acc(_, k, Bool.not(Bool.and(Nat.is_lt(lo, n), Nat.is_lt(lo, t)))), b2n(Nat.is_lt(Nat.add(lo, C.shift(w, _)), Nat.add(t, C.shift(w, k))))) == b2n(Nat.is_lt(Nat.add(lo, C.shift(w, _)), C.shift(w, 1n+k))) : Nat}
114	      %Equal.sym(Bool, Bool.and(Nat.is_lt(lo, n), Nat.is_lt(lo, t)), Nat.is_lt(lo, t), and_lt(lo, n, t, ht, Nat.is_lt(lo, t), {==})) : {Nat.add(hit_acc(k, k, Bool.not(_)), b2n(Nat.is_lt(Nat.add(lo, C.shift(w, k)), Nat.add(t, C.shift(w, k))))) == b2n(Nat.is_lt(Nat.add(lo, C.shift(w, k)), C.shift(w, 1n+k))) : Nat}
115	      term_eq(w, k, t, lo, hlo)
116	    case False{} False{}:
117	      term_gt(w, k, t, lo, hi, Bool.not(Bool.and(Nat.is_lt(lo, n), Nat.is_lt(lo, t))), hab, N.lt_or_eq(k, hi, N.not_lt_le(hi, k, hlt), N.is_eq_sym_false(hi, k, heq)))
118	
119	# THEOREM (one output, Lemire's branch): with t = 2^w mod n < n and
120	# k 2^w + t <= (k + 1) 2^w, y = x n draws k exactly when
121	# k 2^w + t <= y < (k + 1) 2^w
122	def term(+w: Nat, +n: Nat, +k: Nat, +y: Nat, +ht: {Nat.is_lt(SR.pow2mod(w, n), n) == True{} : Bool}, +hab: {Nat.is_le(Nat.add(SR.pow2mod(w, n), C.shift(w, k)), C.shift(w, 1n+k)) == True{} : Bool}) -> {Nat.add(SR.hit(SR.lemire(w, n, C.high(w, y), C.low(w, y)), k), b2n(Nat.is_lt(y, Nat.add(SR.pow2mod(w, n), C.shift(w, k))))) == b2n(Nat.is_lt(y, C.shift(w, 1n+k))) : Nat}:
123	  +lo = C.low(w, y)
124	  +hi = C.high(w, y)
125	  +t = SR.pow2mod(w, n)
126	  %Equal.sym(Nat, y, Nat.add(lo, C.shift(w, hi)), WI.low_high(w, y)) : {Nat.add(SR.hit(SR.lemire(w, n, hi, lo), k), b2n(Nat.is_lt(_, Nat.add(t, C.shift(w, k))))) == b2n(Nat.is_lt(_, C.shift(w, 1n+k))) : Nat}
127	  term_cases(w, n, k, t, lo, hi, WI.low_fits(w, y), ht, hab, Nat.is_lt(hi, k), Nat.is_eq(hi, k), {==}, {==})
128	
129	# ---- summing over the outputs x < N ----
130	
131	# #{x < N : x n < c}
132	def cnt_lt(+n: Nat, +c: Nat, N: Nat) -> Nat:
133	  match N:
134	    case 0n:
135	      0n
136	    case 1n+ +x:
137	      Nat.add(cnt_lt(n, c, x), b2n(Nat.is_lt(Nat.mul(x, n), c)))
138	
139	# the draw is Lemire's when n is not a power of two
140	def draw_lemire(+w: Nat, +m: Nat, +x: Nat, +hp: {Nat.is_eq(SR.and_bits(w, 1n+m, m), 0n) == False{} : Bool}) -> {SR.draw(w, x, 1n+m) == SR.lemire(w, 1n+m, C.high(w, Nat.mul(x, 1n+m)), C.low(w, Nat.mul(x, 1n+m))) : Maybe<&2, Nat>}:
141	  %Equal.sym(Bool, Nat.is_eq(SR.and_bits(w, 1n+m, m), 0n), False{}, hp) : {SR.draw_pos(w, x, 1n+m, m, _) == SR.lemire(w, 1n+m, C.high(w, Nat.mul(x, 1n+m)), C.low(w, Nat.mul(x, 1n+m))) : Maybe<&2, Nat>}
142	  {==}
143	
144	# (a + h) + (b + e) == (a + b) + (h + e)
145	def add4(+a: Nat, +h: Nat, +b: Nat, +e: Nat) -> {Nat.add(Nat.add(a, h), Nat.add(b, e)) == Nat.add(Nat.add(a, b), Nat.add(h, e)) : Nat}:
146	  Equal.trans(Nat, Nat.add(Nat.add(a, h), Nat.add(b, e)), Nat.add(a, Nat.add(h, Nat.add(b, e))), Nat.add(Nat.add(a, b), Nat.add(h, e)),
147	    N.add_assoc(a, h, Nat.add(b, e)),
148	    Equal.trans(Nat, Nat.add(a, Nat.add(h, Nat.add(b, e))), Nat.add(a, Nat.add(b, Nat.add(h, e))), Nat.add(Nat.add(a, b), Nat.add(h, e)),
149	      Equal.cong(Nat, Nat, z => Nat.add(a, z), Nat.add(h, Nat.add(b, e)), Nat.add(b, Nat.add(h, e)),
150	        Equal.trans(Nat, Nat.add(h, Nat.add(b, e)), Nat.add(Nat.add(h, b), e), Nat.add(b, Nat.add(h, e)),
151	          Equal.sym(Nat, Nat.add(Nat.add(h, b), e), Nat.add(h, Nat.add(b, e)), N.add_assoc(h, b, e)),
152	          Equal.trans(Nat, Nat.add(Nat.add(h, b), e), Nat.add(Nat.add(b, h), e), Nat.add(b, Nat.add(h, e)),
153	            Equal.cong(Nat, Nat, z => Nat.add(z, e), Nat.add(h, b), Nat.add(b, h), N.add_comm(h, b)),
154	            N.add_assoc(b, h, e)))),
155	      Equal.sym(Nat, Nat.add(Nat.add(a, b), Nat.add(h, e)), Nat.add(a, Nat.add(b, Nat.add(h, e))), N.add_assoc(a, b, Nat.add(h, e)))))
156	
157	# the Lemire outputs below N plus those with x n below the lower threshold
158	# are those with x n below the upper one
159	def sum_lemire(+w: Nat, +m: Nat, +k: Nat, +N: Nat, +hp: {Nat.is_eq(SR.and_bits(w, 1n+m, m), 0n) == False{} : Bool}, +ht: {Nat.is_lt(SR.pow2mod(w, 1n+m), 1n+m) == True{} : Bool}, +hab: {Nat.is_le(Nat.add(SR.pow2mod(w, 1n+m), C.shift(w, k)), C.shift(w, 1n+k)) == True{} : Bool}) -> {Nat.add(SR.count(w, 1n+m, k, N), cnt_lt(1n+m, Nat.add(SR.pow2mod(w, 1n+m), C.shift(w, k)), N)) == cnt_lt(1n+m, C.shift(w, 1n+k), N) : Nat}:
160	  match N:
161	    case 0n:
162	      {==}
163	    case 1n+ +x:
164	      +a = Nat.add(SR.pow2mod(w, 1n+m), C.shift(w, k))
165	      +b = C.shift(w, 1n+k)
166	      +h = SR.hit(SR.draw(w, x, 1n+m), k)
167	      +ea = b2n(Nat.is_lt(Nat.mul(x, 1n+m), a))
168	      Equal.trans(Nat, Nat.add(Nat.add(SR.count(w, 1n+m, k, x), h), Nat.add(cnt_lt(1n+m, a, x), ea)), Nat.add(Nat.add(SR.count(w, 1n+m, k, x), cnt_lt(1n+m, a, x)), Nat.add(h, ea)), Nat.add(cnt_lt(1n+m, b, x), b2n(Nat.is_lt(Nat.mul(x, 1n+m), b))),
169	        add4(SR.count(w, 1n+m, k, x), h, cnt_lt(1n+m, a, x), ea),
170	        Equal.trans(Nat, Nat.add(Nat.add(SR.count(w, 1n+m, k, x), cnt_lt(1n+m, a, x)), Nat.add(h, ea)), Nat.add(cnt_lt(1n+m, b, x), Nat.add(h, ea)), Nat.add(cnt_lt(1n+m, b, x), b2n(Nat.is_lt(Nat.mul(x, 1n+m), b))),
171	          Equal.cong(Nat, Nat, z => Nat.add(z, Nat.add(h, ea)), Nat.add(SR.count(w, 1n+m, k, x), cnt_lt(1n+m, a, x)), cnt_lt(1n+m, b, x), sum_lemire(w, m, k, x, hp, ht, hab)),
172	          Equal.cong(Nat, Nat, z => Nat.add(cnt_lt(1n+m, b, x), z), Nat.add(h, ea), b2n(Nat.is_lt(Nat.mul(x, 1n+m), b)),
173	            %Equal.sym(Maybe<&2, Nat>, SR.draw(w, x, 1n+m), SR.lemire(w, 1n+m, C.high(w, Nat.mul(x, 1n+m)), C.low(w, Nat.mul(x, 1n+m))), draw_lemire(w, m, x, hp)) : {Nat.add(SR.hit(_, k), ea) == b2n(Nat.is_lt(Nat.mul(x, 1n+m), b)) : Nat}
174	            term(w, 1n+m, k, Nat.mul(x, 1n+m), ht, hab))))
175	
176	# ---- counting x n < c: min(N, ceil(c / n)) ----
177	
178	# a + 1 <= b is a < b
179	def le_succ_lt(+a: Nat, +b: Nat) -> {Nat.is_le(1n+a, b) == Nat.is_lt(a, b) : Bool}:
180	  match a b:
181	    case 0n 0n:
182	      {==}
183	    case 0n 1n+q:
184	      N.zero_le(q)
185	    case 1n+p 0n:
186	      {==}
187	    case 1n+p 1n+q:
188	      le_succ_lt(p, q)
189	
190	# x n < c is x < ceil(c / n) = (c + n - 1) / n, n = 1 + bp
191	def lt_ceil(+bp: Nat, +x: Nat, +c: Nat) -> {Nat.is_lt(Nat.mul(x, 1n+bp), c) == Nat.is_lt(x, Nat.div(Nat.add(c, bp), 1n+bp)) : Bool}:
192	  +X = Nat.mul(x, 1n+bp)
193	  +D = Nat.div(Nat.add(c, bp), 1n+bp)
194	  Equal.trans(Bool, Nat.is_lt(X, c), Nat.is_le(1n+X, c), Nat.is_lt(x, D),
195	    Equal.sym(Bool, Nat.is_le(1n+X, c), Nat.is_lt(X, c), le_succ_lt(X, c)),
196	    Equal.trans(Bool, Nat.is_le(1n+X, c), Nat.is_le(Nat.add(bp, 1n+X), Nat.add(bp, c)), Nat.is_lt(x, D),
197	      Equal.sym(Bool, Nat.is_le(Nat.add(bp, 1n+X), Nat.add(bp, c)), Nat.is_le(1n+X, c), AR.le_add_cancel(bp, 1n+X, c)),
198	      Equal.trans(Bool, Nat.is_le(Nat.add(bp, 1n+X), Nat.add(bp, c)), Nat.is_le(1n+Nat.add(bp, X), Nat.add(c, bp)), Nat.is_lt(x, D),
199	        %N.add_succ(bp, X) : {Nat.is_le(Nat.add(bp, 1n+X), Nat.add(bp, c)) == Nat.is_le(_, Nat.add(c, bp)) : Bool}
200	        %N.add_comm(bp, c) : {Nat.is_le(Nat.add(bp, 1n+X), Nat.add(bp, c)) == Nat.is_le(Nat.add(bp, 1n+X), _) : Bool}
201	        {==},
202	        Equal.trans(Bool, Nat.is_le(1n+Nat.add(bp, X), Nat.add(c, bp)), Nat.is_le(1n+x, D), Nat.is_lt(x, D),
203	          Equal.sym(Bool, Nat.is_le(1n+x, D), Nat.is_le(1n+Nat.add(bp, X), Nat.add(c, bp)), AR.le_div(bp, 1n+x, Nat.add(c, bp))),
204	          le_succ_lt(x, D)))))
205	
206	def min_zero(+n: Nat) -> {Nat.min(n, 0n) == 0n : Nat}:
207	  match n:
208	    case 0n:
209	      {==}
210	    case 1n+p:
211	      {==}
212	
213	# min(N + 1, C) = min(N, C) + [N < C]
214	def min_succ(+n: Nat, +c: Nat) -> {Nat.min(1n+n, c) == Nat.add(Nat.min(n, c), b2n(Nat.is_lt(n, c))) : Nat}:
215	  match n c:
216	    case _ 0n:
217	      %Equal.sym(Nat, Nat.min(n, 0n), 0n, min_zero(n)) : {0n == Nat.add(_, b2n(Nat.is_lt(n, 0n))) : Nat}
218	      %Equal.sym(Bool, Nat.is_lt(n, 0n), False{}, N.not_lt_zero(n)) : {0n == Nat.add(0n, b2n(_)) : Nat}
219	      {==}
220	    case 0n 1n+q:
221	      {==}
222	    case 1n+p 1n+q:
223	      Equal.cong(Nat, Nat, z => 1n+z, Nat.min(1n+p, q), Nat.add(Nat.min(p, q), b2n(Nat.is_lt(p, q))), min_succ(p, q))
224	
225	# THEOREM: #{x < N : x n < c} = min(N, ceil(c / n))
226	def cnt_min(+bp: Nat, +c: Nat, +N: Nat) -> {cnt_lt(1n+bp, c, N) == Nat.min(N, Nat.div(Nat.add(c, bp), 1n+bp)) : Nat}:
227	  match N:
228	    case 0n:
229	      {==}
230	    case 1n+ +x:
231	      +D = Nat.div(Nat.add(c, bp), 1n+bp)
232	      Equal.trans(Nat, Nat.add(cnt_lt(1n+bp, c, x), b2n(Nat.is_lt(Nat.mul(x, 1n+bp), c))), Nat.add(Nat.min(x, D), b2n(Nat.is_lt(x, D))), Nat.min(1n+x, D),
233	        %Equal.sym(Nat, cnt_lt(1n+bp, c, x), Nat.min(x, D), cnt_min(bp, c, x)) : {Nat.add(_, b2n(Nat.is_lt(Nat.mul(x, 1n+bp), c))) == Nat.add(Nat.min(x, D), b2n(Nat.is_lt(x, D))) : Nat}
234	        %Equal.sym(Bool, Nat.is_lt(Nat.mul(x, 1n+bp), c), Nat.is_lt(x, D), lt_ceil(bp, x, c)) : {Nat.add(Nat.min(x, D), b2n(_)) == Nat.add(Nat.min(x, D), b2n(Nat.is_lt(x, D))) : Nat}
235	        {==},
236	        Equal.sym(Nat, Nat.min(1n+x, D), Nat.add(Nat.min(x, D), b2n(Nat.is_lt(x, D))), min_succ(x, D)))
237	
238	# ---- Lemire's branch: the count ----
239	
240	def dbl_add(+z: Nat) -> {Nat.double(z) == Nat.add(z, z) : Nat}:
241	  match z:
242	    case 0n:
243	      {==}
244	    case 1n+p:
245	      %Equal.sym(Nat, Nat.add(p, 1n+p), 1n+Nat.add(p, p), N.add_succ(p, p)) : {2n+Nat.double(p) == 1n+_ : Nat}
246	      Equal.cong(Nat, Nat, x => 2n+x, Nat.double(p), Nat.add(p, p), dbl_add(p))
247	
248	# Go's -n % n: 2^w mod n
249	def pm(+w: Nat, +m: Nat) -> {SR.pow2mod(w, 1n+m) == Nat.mod(C.shift(w, 1n), 1n+m) : Nat}:
250	  match w:
251	    case 0n:
252	      {==}
253	    case 1n+k:
254	      +S = C.shift(k, 1n)
255	      %Equal.sym(Nat, SR.pow2mod(k, 1n+m), Nat.mod(S, 1n+m), pm(k, m)) : {Nat.mod(Nat.double(_), 1n+m) == Nat.mod(Nat.double(S), 1n+m) : Nat}
256	      %Equal.sym(Nat, Nat.double(Nat.mod(S, 1n+m)), Nat.add(Nat.mod(S, 1n+m), Nat.mod(S, 1n+m)), dbl_add(Nat.mod(S, 1n+m))) : {Nat.mod(_, 1n+m) == Nat.mod(Nat.double(S), 1n+m) : Nat}
257	      %Equal.sym(Nat, Nat.double(S), Nat.add(S, S), dbl_add(S)) : {Nat.mod(Nat.add(Nat.mod(S, 1n+m), Nat.mod(S, 1n+m)), 1n+m) == Nat.mod(_, 1n+m) : Nat}
258	      Equal.trans(Nat, Nat.mod(Nat.add(Nat.mod(S, 1n+m), Nat.mod(S, 1n+m)), 1n+m), Nat.mod(Nat.add(S, Nat.mod(S, 1n+m)), 1n+m), Nat.mod(Nat.add(S, S), 1n+m),
259	        AR.mod_add_l(m, S, Nat.mod(S, 1n+m)),
260	        AR.mod_add_r(m, S, S))
261	
262	# (q n + r) / n = q + r / n
263	def div_add_mul(+q: Nat, +m: Nat, +r: Nat) -> {Nat.div(Nat.add(Nat.mul(q, 1n+m), r), 1n+m) == Nat.add(q, Nat.div(r, 1n+m)) : Nat}:
264	  +rq = Nat.div(r, 1n+m)
265	  +rr = Nat.mod(r, 1n+m)
266	  %Equal.sym(Nat, r, Nat.add(Nat.mul(rq, 1n+m), rr), AR.dm_eq(m, r)) : {Nat.div(Nat.add(Nat.mul(q, 1n+m), _), 1n+m) == Nat.add(q, rq) : Nat}
267	  %N.add_assoc(Nat.mul(q, 1n+m), Nat.mul(rq, 1n+m), rr) : {Nat.div(_, 1n+m) == Nat.add(q, rq) : Nat}
268	  %NA.mul_add_right(q, rq, 1n+m) : {Nat.div(Nat.add(_, rr), 1n+m) == Nat.add(q, rq) : Nat}
269	  AR.div_of(Nat.add(q, rq), m, rr, AR.dm_lt(m, r))
270	
271	def min_le(+n: Nat, +c: Nat, +h: {Nat.is_le(c, n) == True{} : Bool}) -> {Nat.min(n, c) == c : Nat}:
272	  match n c:
273	    case _ 0n:
274	      min_zero(n)
275	    case 0n 1n+q:
276	      Empty.absurd({Nat.min(0n, 1n+q) == 1n+q : Nat}, L.false_true(h))
277	    case 1n+p 1n+q:
278	      Equal.cong(Nat, Nat, z => 1n+z, Nat.min(p, q), q, min_le(p, q, h))
279	
280	# a <= b and c < d give a + c < b + d
281	def lt_add2(+a: Nat, +b: Nat, +c: Nat, +d: Nat, +hab: {Nat.is_le(a, b) == True{} : Bool}, +hcd: {Nat.is_lt(c, d) == True{} : Bool}) -> {Nat.is_lt(Nat.add(a, c), Nat.add(b, d)) == True{} : Bool}:
282	  N.le_lt_trans(Nat.add(a, c), Nat.add(b, c), Nat.add(b, d), WI.le_add_r(a, b, c, hab), N.lt_add_left(c, d, b, hcd))
283	
284	# ceil((k + 1) 2^w / n) <= 2^w for k < n
285	def cb_le(+w: Nat, +m: Nat, +k: Nat, +hk: {Nat.is_lt(k, 1n+m) == True{} : Bool}) -> {Nat.is_le(Nat.div(Nat.add(C.shift(w, 1n+k), m), 1n+m), C.shift(w, 1n)) == True{} : Bool}:
286	  +W = C.shift(w, 1n)
287	  +b = C.shift(w, 1n+k)
288	  +hb = Equal.trans(Bool, Nat.is_le(b, Nat.mul(W, 1n+m)), Nat.is_le(Nat.mul(1n+k, W), Nat.mul(1n+m, W)), True{},
289	    %WI.shift_mul(w, 1n+k) : {Nat.is_le(b, Nat.mul(W, 1n+m)) == Nat.is_le(_, Nat.mul(1n+m, W)) : Bool}
290	    %NA.mul_comm(W, 1n+m) : {Nat.is_le(b, Nat.mul(W, 1n+m)) == Nat.is_le(b, _) : Bool}
291	    {==},
292	    LA.mul_le(1n+k, 1n+m, W, N.lt_succ_le_succ(k, 1n+m, hk)))
293	  +hlt = lt_add2(b, Nat.mul(W, 1n+m), m, 1n+m, hb, N.lt_succ(m))
294	  +hlt2 = Equal.trans(Bool, Nat.is_lt(Nat.add(b, m), Nat.mul(1n+W, 1n+m)), Nat.is_lt(Nat.add(b, m), Nat.add(Nat.mul(W, 1n+m), 1n+m)), True{},
295	    %N.add_comm(1n+m, Nat.mul(W, 1n+m)) : {Nat.is_lt(Nat.add(b, m), Nat.mul(1n+W, 1n+m)) == Nat.is_lt(Nat.add(b, m), _) : Bool}
296	    {==},
297	    hlt)
298	  +hnot = Equal.trans(Bool, Nat.is_le(1n+W, Nat.div(Nat.add(b, m), 1n+m)), Nat.is_le(Nat.mul(1n+W, 1n+m), Nat.add(b, m)), False{},
299	    AR.le_div(m, 1n+W, Nat.add(b, m)),
300	    N.lt_not_le(Nat.add(b, m), Nat.mul(1n+W, 1n+m), hlt2))
301	  N.lt_succ_le(Nat.div(Nat.add(b, m), 1n+m), W, N.not_le_lt(1n+W, Nat.div(Nat.add(b, m), 1n+m), hnot))
302	
303	# t <= 2^w for t < n <= 2^w
304	def hab_of(+w: Nat, +m: Nat, +k: Nat, +hw: {C.fits(w, 1n+m) == True{} : Bool}) -> {Nat.is_le(Nat.add(SR.pow2mod(w, 1n+m), C.shift(w, k)), C.shift(w, 1n+k)) == True{} : Bool}:
305	  +W = C.shift(w, 1n)
306	  +t = SR.pow2mod(w, 1n+m)
307	  +ht = N.lt_le(t, W, N.lt_trans(t, 1n+m, W,
308	    %Equal.sym(Nat, SR.pow2mod(w, 1n+m), Nat.mod(W, 1n+m), pm(w, m)) : {Nat.is_lt(_, 1n+m) == True{} : Bool}
309	    AR.dm_lt(m, W),
310	    WI.lt_one(w, 1n, {==}, 1n+m, hw)))
311	  %Equal.sym(Nat, C.shift(w, 1n+k), Nat.add(W, C.shift(w, k)), WI.shift_add(w, 1n, k)) : {Nat.is_le(Nat.add(t, C.shift(w, k)), _) == True{} : Bool}
312	  WI.le_add_r(t, W, C.shift(w, k), ht)
313	
314	# (k + 1) 2^w + m = q n + (t + k 2^w + m), q n + t = 2^w
315	def bm_eq(+w: Nat, +m: Nat, +k: Nat) -> {Nat.add(C.shift(w, 1n+k), m) == Nat.add(Nat.mul(Nat.div(C.shift(w, 1n), 1n+m), 1n+m), Nat.add(Nat.add(SR.pow2mod(w, 1n+m), C.shift(w, k)), m)) : Nat}:
316	  +W = C.shift(w, 1n)
317	  +S = C.shift(w, k)
318	  +qn = Nat.mul(Nat.div(W, 1n+m), 1n+m)
319	  +t = SR.pow2mod(w, 1n+m)
320	  %Equal.sym(Nat, C.shift(w, 1n+k), Nat.add(W, S), WI.shift_add(w, 1n, k)) : {Nat.add(_, m) == Nat.add(qn, Nat.add(Nat.add(t, S), m)) : Nat}
321	  %Equal.sym(Nat, W, Nat.add(qn, Nat.mod(W, 1n+m)), AR.dm_eq(m, W)) : {Nat.add(Nat.add(_, S), m) == Nat.add(qn, Nat.add(Nat.add(t, S), m)) : Nat}
322	  %pm(w, m) : {Nat.add(Nat.add(Nat.add(qn, _), S), m) == Nat.add(qn, Nat.add(Nat.add(t, S), m)) : Nat}
323	  %Equal.sym(Nat, Nat.add(Nat.add(qn, t), S), Nat.add(qn, Nat.add(t, S)), N.add_assoc(qn, t, S)) : {Nat.add(_, m) == Nat.add(qn, Nat.add(Nat.add(t, S), m)) : Nat}
324	  N.add_assoc(qn, Nat.add(t, S), m)
325	
326	# THEOREM (Lemire's branch): floor(2^w / n) outputs draw k
327	def lemire_count(+w: Nat, +m: Nat, +k: Nat, +hp: {Nat.is_eq(SR.and_bits(w, 1n+m, m), 0n) == False{} : Bool}, +hw: {C.fits(w, 1n+m) == True{} : Bool}, +hk: {Nat.is_lt(k, 1n+m) == True{} : Bool}) -> {SR.count(w, 1n+m, k, C.shift(w, 1n)) == Nat.div(C.shift(w, 1n), 1n+m) : Nat}:
328	  +W = C.shift(w, 1n)
329	  +t = SR.pow2mod(w, 1n+m)
330	  +a = Nat.add(t, C.shift(w, k))
331	  +b = C.shift(w, 1n+k)
332	  +q = Nat.div(W, 1n+m)
333	  +Ca = Nat.div(Nat.add(a, m), 1n+m)
334	  +Cb = Nat.div(Nat.add(b, m), 1n+m)
335	  +cnt = SR.count(w, 1n+m, k, W)
336	  +ht = Equal.trans(Bool, Nat.is_lt(t, 1n+m), Nat.is_lt(Nat.mod(W, 1n+m), 1n+m), True{}, Equal.cong(Nat, Bool, z => Nat.is_lt(z, 1n+m), t, Nat.mod(W, 1n+m), pm(w, m)), AR.dm_lt(m, W))
337	  +hab = hab_of(w, m, k, hw)
338	  +hcb = cb_le(w, m, k, hk)
339	  +bq = Equal.trans(Nat, Cb, Nat.div(Nat.add(Nat.mul(q, 1n+m), Nat.add(a, m)), 1n+m), Nat.add(q, Ca),
340	    Equal.cong(Nat, Nat, z => Nat.div(z, 1n+m), Nat.add(b, m), Nat.add(Nat.mul(q, 1n+m), Nat.add(a, m)), bm_eq(w, m, k)),
341	    div_add_mul(q, m, Nat.add(a, m)))
342	  +hca = Equal.trans(Bool, Nat.is_le(Ca, Cb), Nat.is_le(Ca, Nat.add(q, Ca)), True{}, Equal.cong(Nat, Bool, z => Nat.is_le(Ca, z), Cb, Nat.add(q, Ca), bq), le_plus(q, Ca))
343	  +ma = min_le(W, Ca, N.le_trans(Ca, Cb, W, hca, hcb))
344	  +mb = min_le(W, Cb, hcb)
345	  UA.add_cancel_r(cnt, q, Ca,
346	    Equal.trans(Nat, Nat.add(cnt, Ca), Nat.add(cnt, cnt_lt(1n+m, a, W)), Nat.add(q, Ca),
347	      Equal.cong(Nat, Nat, z => Nat.add(cnt, z), Ca, cnt_lt(1n+m, a, W),
348	        Equal.trans(Nat, Ca, Nat.min(W, Ca), cnt_lt(1n+m, a, W), Equal.sym(Nat, Nat.min(W, Ca), Ca, ma), Equal.sym(Nat, cnt_lt(1n+m, a, W), Nat.min(W, Ca), cnt_min(m, a, W)))),
349	      Equal.trans(Nat, Nat.add(cnt, cnt_lt(1n+m, a, W)), cnt_lt(1n+m, b, W), Nat.add(q, Ca),
350	        sum_lemire(w, m, k, W, hp, ht, hab),
351	        Equal.trans(Nat, cnt_lt(1n+m, b, W), Nat.min(W, Cb), Nat.add(q, Ca),
352	          cnt_min(m, b, W),
353	          Equal.trans(Nat, Nat.min(W, Cb), Cb, Nat.add(q, Ca), mb, bq)))))
354	
355	# ---- the mask branch: n = 2^j, x & (n - 1) ----
356	
357	# #{x < N : and_bits(w, x, m) == k}
358	def cm(+w: Nat, +m: Nat, +k: Nat, N: Nat) -> Nat:
359	  match N:
360	    case 0n:
361	      0n
362	    case 1n+ +x:
363	      Nat.add(cm(w, m, k, x), b2n(Nat.is_eq(SR.and_bits(w, x, m), k)))
364	
365	# the same count taken over the pairs 2y, 2y + 1 with y < N
366	def cm2(+w: Nat, +m: Nat, +k: Nat, N: Nat) -> Nat:
367	  match N:
368	    case 0n:
369	      0n
370	    case 1n+ +y:
371	      Nat.add(cm2(w, m, k, y), Nat.add(b2n(Nat.is_eq(SR.and_bits(w, Nat.double(y), m), k)), b2n(Nat.is_eq(SR.and_bits(w, 1n+Nat.double(y), m), k))))
372	
373	def split(+w: Nat, +m: Nat, +k: Nat, +N: Nat) -> {cm(w, m, k, Nat.double(N)) == cm2(w, m, k, N) : Nat}:
374	  match N:
375	    case 0n:
376	      {==}
377	    case 1n+ +y:
378	      +f0 = b2n(Nat.is_eq(SR.and_bits(w, Nat.double(y), m), k))
379	      +f1 = b2n(Nat.is_eq(SR.and_bits(w, 1n+Nat.double(y), m), k))
380	      Equal.trans(Nat, Nat.add(Nat.add(cm(w, m, k, Nat.double(y)), f0), f1), Nat.add(cm(w, m, k, Nat.double(y)), Nat.add(f0, f1)), Nat.add(cm2(w, m, k, y), Nat.add(f0, f1)),
381	        N.add_assoc(cm(w, m, k, Nat.double(y)), f0, f1),
382	        Equal.cong(Nat, Nat, z => Nat.add(z, Nat.add(f0, f1)), cm(w, m, k, Nat.double(y)), cm2(w, m, k, y), split(w, m, k, y)))
383	
384	# [2a == k] + [2a + 1 == k] == [a == k / 2]
385	def pair_eq(+a: Nat, +k: Nat) -> {Nat.add(b2n(Nat.is_eq(Nat.double(a), k)), b2n(Nat.is_eq(1n+Nat.double(a), k))) == b2n(Nat.is_eq(a, C.half(k))) : Nat}:
386	  match a k:
387	    case 0n 0n:
388	      {==}
389	    case 0n 1n:
390	      {==}
391	    case 0n 2n+q:
392	      {==}
393	    case 1n+p 0n:
394	      {==}
395	    case 1n+p 1n:
396	      {==}
397	    case 1n+p 2n+q:
398	      pair_eq(p, q)
399	
400	def bit_dbl0(+y: Nat) -> {C.bit(Nat.double(y)) == 0n : Nat}:
401	  WI.bit_dbl(0n, y)
402	
403	def half_dbl0(+y: Nat) -> {C.half(Nat.double(y)) == y : Nat}:
404	  WI.half_dbl(0n, y)
405	
406	def bit_dbl1(+y: Nat) -> {C.bit(1n+Nat.double(y)) == 1n : Nat}:
407	  WI.bit_dbl(1n, y)
408	
409	def half_dbl1(+y: Nat) -> {C.half(1n+Nat.double(y)) == y : Nat}:
410	  WI.half_dbl(1n, y)
411	
412	# with m = 2h + 1: the pair 2y, 2y + 1 draws k once exactly when y draws k / 2
413	# one width lower with h
414	def pair_and(+v: Nat, +h: Nat, +k: Nat, +y: Nat) -> {Nat.add(b2n(Nat.is_eq(SR.and_bits(1n+v, Nat.double(y), 1n+Nat.double(h)), k)), b2n(Nat.is_eq(SR.and_bits(1n+v, 1n+Nat.double(y), 1n+Nat.double(h)), k))) == b2n(Nat.is_eq(SR.and_bits(v, y, h), C.half(k))) : Nat}:
415	  +M = {1n+Nat.double(h) : Nat}
416	  %Equal.sym(Nat, C.bit(M), 1n, bit_dbl1(h)) : {Nat.add(b2n(Nat.is_eq(Nat.add(Nat.mul(C.bit(Nat.double(y)), _), Nat.double(SR.and_bits(v, C.half(Nat.double(y)), C.half(M)))), k)), b2n(Nat.is_eq(Nat.add(Nat.mul(C.bit(1n+Nat.double(y)), _), Nat.double(SR.and_bits(v, C.half(1n+Nat.double(y)), C.half(M)))), k))) == b2n(Nat.is_eq(SR.and_bits(v, y, h), C.half(k))) : Nat}
417	  %Equal.sym(Nat, C.half(M), h, half_dbl1(h)) : {Nat.add(b2n(Nat.is_eq(Nat.add(Nat.mul(C.bit(Nat.double(y)), 1n), Nat.double(SR.and_bits(v, C.half(Nat.double(y)), _))), k)), b2n(Nat.is_eq(Nat.add(Nat.mul(C.bit(1n+Nat.double(y)), 1n), Nat.double(SR.and_bits(v, C.half(1n+Nat.double(y)), _))), k))) == b2n(Nat.is_eq(SR.and_bits(v, y, h), C.half(k))) : Nat}
418	  %Equal.sym(Nat, C.bit(Nat.double(y)), 0n, bit_dbl0(y)) : {Nat.add(b2n(Nat.is_eq(Nat.add(Nat.mul(_, 1n), Nat.double(SR.and_bits(v, C.half(Nat.double(y)), h))), k)), b2n(Nat.is_eq(Nat.add(Nat.mul(C.bit(1n+Nat.double(y)), 1n), Nat.double(SR.and_bits(v, C.half(1n+Nat.double(y)), h))), k))) == b2n(Nat.is_eq(SR.and_bits(v, y, h), C.half(k))) : Nat}
419	  %Equal.sym(Nat, C.half(Nat.double(y)), y, half_dbl0(y)) : {Nat.add(b2n(Nat.is_eq(Nat.add(Nat.mul(0n, 1n), Nat.double(SR.and_bits(v, _, h))), k)), b2n(Nat.is_eq(Nat.add(Nat.mul(C.bit(1n+Nat.double(y)), 1n), Nat.double(SR.and_bits(v, C.half(1n+Nat.double(y)), h))), k))) == b2n(Nat.is_eq(SR.and_bits(v, y, h), C.half(k))) : Nat}
420	  %Equal.sym(Nat, C.bit(1n+Nat.double(y)), 1n, bit_dbl1(y)) : {Nat.add(b2n(Nat.is_eq(Nat.double(SR.and_bits(v, y, h)), k)), b2n(Nat.is_eq(Nat.add(Nat.mul(_, 1n), Nat.double(SR.and_bits(v, C.half(1n+Nat.double(y)), h))), k))) == b2n(Nat.is_eq(SR.and_bits(v, y, h), C.half(k))) : Nat}
421	  %Equal.sym(Nat, C.half(1n+Nat.double(y)), y, half_dbl1(y)) : {Nat.add(b2n(Nat.is_eq(Nat.double(SR.and_bits(v, y, h)), k)), b2n(Nat.is_eq(Nat.add(1n, Nat.double(SR.and_bits(v, _, h))), k))) == b2n(Nat.is_eq(SR.and_bits(v, y, h), C.half(k))) : Nat}
422	  pair_eq(SR.and_bits(v, y, h), k)
423	
424	# the pairs one width up count what the lower width counts
425	def pairs(+v: Nat, +h: Nat, +k: Nat, +N: Nat) -> {cm2(1n+v, 1n+Nat.double(h), k, N) == cm(v, h, C.half(k), N) : Nat}:
426	  match N:
427	    case 0n:
428	      {==}
429	    case 1n+ +y:
430	      Equal.trans(Nat, Nat.add(cm2(1n+v, 1n+Nat.double(h), k, y), Nat.add(b2n(Nat.is_eq(SR.and_bits(1n+v, Nat.double(y), 1n+Nat.double(h)), k)), b2n(Nat.is_eq(SR.and_bits(1n+v, 1n+Nat.double(y), 1n+Nat.double(h)), k)))), Nat.add(cm(v, h, C.half(k), y), Nat.add(b2n(Nat.is_eq(SR.and_bits(1n+v, Nat.double(y), 1n+Nat.double(h)), k)), b2n(Nat.is_eq(SR.and_bits(1n+v, 1n+Nat.double(y), 1n+Nat.double(h)), k)))), Nat.add(cm(v, h, C.half(k), y), b2n(Nat.is_eq(SR.and_bits(v, y, h), C.half(k)))),
431	        Equal.cong(Nat, Nat, z => Nat.add(z, Nat.add(b2n(Nat.is_eq(SR.and_bits(1n+v, Nat.double(y), 1n+Nat.double(h)), k)), b2n(Nat.is_eq(SR.and_bits(1n+v, 1n+Nat.double(y), 1n+Nat.double(h)), k)))), cm2(1n+v, 1n+Nat.double(h), k, y), cm(v, h, C.half(k), y), pairs(v, h, k, y)),
432	        Equal.cong(Nat, Nat, z => Nat.add(cm(v, h, C.half(k), y), z), Nat.add(b2n(Nat.is_eq(SR.and_bits(1n+v, Nat.double(y), 1n+Nat.double(h)), k)), b2n(Nat.is_eq(SR.and_bits(1n+v, 1n+Nat.double(y), 1n+Nat.double(h)), k))), b2n(Nat.is_eq(SR.and_bits(v, y, h), C.half(k))), pair_and(v, h, k, y)))
433	
434	# x & 0 == 0
435	def and0(+w: Nat, +x: Nat) -> {SR.and_bits(w, x, 0n) == 0n : Nat}:
436	  match w:
437	    case 0n:
438	      {==}
439	    case 1n+v:
440	      %Equal.sym(Nat, Nat.mul(C.bit(x), 0n), 0n, NA.mul_zero(C.bit(x))) : {Nat.add(_, Nat.double(SR.and_bits(v, C.half(x), 0n))) == 0n : Nat}
441	      %Equal.sym(Nat, SR.and_bits(v, C.half(x), 0n), 0n, and0(v, C.half(x))) : {Nat.add(0n, Nat.double(_)) == 0n : Nat}
442	      {==}
443	
444	# n = 1: every output draws 0
445	def cm0(+w: Nat, +N: Nat) -> {cm(w, 0n, 0n, N) == N : Nat}:
446	  match N:
447	    case 0n:
448	      {==}
449	    case 1n+ +x:
450	      %Equal.sym(Nat, SR.and_bits(w, x, 0n), 0n, and0(w, x)) : {Nat.add(cm(w, 0n, 0n, x), b2n(Nat.is_eq(_, 0n))) == 1n+x : Nat}
451	      %Equal.sym(Nat, cm(w, 0n, 0n, x), x, cm0(w, x)) : {Nat.add(_, 1n) == 1n+x : Nat}
452	      Equal.trans(Nat, Nat.add(x, 1n), 1n+Nat.add(x, 0n), 1n+x, N.add_succ(x, 0n), Equal.cong(Nat, Nat, z => 1n+z, Nat.add(x, 0n), x, N.add_zero(x)))
453	
454	def cm0k(+w: Nat, +k: Nat, +N: Nat, +hk: {Nat.is_lt(k, 1n) == True{} : Bool}) -> {cm(w, 0n, k, N) == N : Nat}:
455	  match k:
456	    case 0n:
457	      cm0(w, N)
458	    case 1n+q:
459	      Empty.absurd({cm(w, 0n, 1n+q, N) == N : Nat}, N.lt_zero_absurd(q, hk))
460	
461	def div_one(+N: Nat) -> {Nat.div(N, 1n) == N : Nat}:
462	  %Equal.sym(Nat, N, Nat.add(Nat.mul(N, 1n), 0n), Equal.trans(Nat, N, Nat.mul(N, 1n), Nat.add(Nat.mul(N, 1n), 0n), Equal.sym(Nat, Nat.mul(N, 1n), N, NA.mul_one(N)), Equal.sym(Nat, Nat.add(Nat.mul(N, 1n), 0n), Nat.mul(N, 1n), N.add_zero(Nat.mul(N, 1n))))) : {Nat.div(_, 1n) == N : Nat}
463	  AR.div_of(N, 0n, 0n, {==})
464	
465	# 2a / 2n = a / n
466	def div_double(+a: Nat, +h: Nat) -> {Nat.div(Nat.double(a), 2n+Nat.double(h)) == Nat.div(a, 1n+h) : Nat}:
467	  +q = Nat.div(a, 1n+h)
468	  +r = Nat.mod(a, 1n+h)
469	  %Equal.sym(Nat, a, Nat.add(Nat.mul(q, 1n+h), r), AR.dm_eq(h, a)) : {Nat.div(Nat.double(_), 2n+Nat.double(h)) == q : Nat}
470	  %N.add_double(Nat.mul(q, 1n+h), r) : {Nat.div(_, 2n+Nat.double(h)) == q : Nat}
471	  %Equal.sym(Nat, Nat.double(Nat.mul(q, 1n+h)), Nat.mul(q, Nat.double(1n+h)), WI.dbl_mul(q, 1n+h)) : {Nat.div(Nat.add(_, Nat.double(r)), 2n+Nat.double(h)) == q : Nat}
472	  AR.div_of(q, 1n+Nat.double(h), Nat.double(r), N.double_lt(r, 1n+h, AR.dm_lt(h, a)))
473	
474	# bit * bit == bit
475	def bitsq_b(+b: Nat, +h: {Nat.is_le(b, 1n) == True{} : Bool}) -> {Nat.mul(b, b) == b : Nat}:
476	  match b:
477	    case 0n:
478	      {==}
479	    case 1n:
480	      {==}
481	    case 2n+q:
482	      Empty.absurd({Nat.mul(2n+q, 2n+q) == 2n+q : Nat}, L.false_true(h))
483	
484	# x & x is the low w bits of x
485	def and_self(+w: Nat, +x: Nat) -> {SR.and_bits(w, x, x) == C.low(w, x) : Nat}:
486	  match w:
487	    case 0n:
488	      {==}
489	    case 1n+v:
490	      %Equal.sym(Nat, Nat.mul(C.bit(x), C.bit(x)), C.bit(x), bitsq_b(C.bit(x), WI.bit_le1(x))) : {Nat.add(_, Nat.double(SR.and_bits(v, C.half(x), C.half(x)))) == Nat.add(C.bit(x), Nat.double(C.low(v, C.half(x)))) : Nat}
491	      Equal.cong(Nat, Nat, z => Nat.add(C.bit(x), Nat.double(z)), SR.and_bits(v, C.half(x), C.half(x)), C.low(v, C.half(x)), and_self(v, C.half(x)))
492	
493	def dbl0(+z: Nat, +h: {Nat.double(z) == 0n : Nat}) -> {z == 0n : Nat}:
494	  match z:
495	    case 0n:
496	      {==}
497	    case 1n+p:
498	      Empty.absurd({1n+p == 0n : Nat}, N.succ_zero(1n+Nat.double(p), h))
499	
500	# m = bit(m) + 2 half(m)
501	def decomp(+m: Nat, +b: Nat, +h: Nat, +hb: {C.bit(m) == b : Nat}, +hh: {C.half(m) == h : Nat}) -> {m == Nat.add(b, Nat.double(h)) : Nat}:
502	  Equal.trans(Nat, m, Nat.add(C.bit(m), Nat.double(C.half(m))), Nat.add(b, Nat.double(h)), WI.hb(m),
503	    Equal.trans(Nat, Nat.add(C.bit(m), Nat.double(C.half(m))), Nat.add(b, Nat.double(C.half(m))), Nat.add(b, Nat.double(h)),
504	      Equal.cong(Nat, Nat, z => Nat.add(z, Nat.double(C.half(m))), C.bit(m), b, hb),
505	      Equal.cong(Nat, Nat, z => Nat.add(b, Nat.double(z)), C.half(m), h, hh)))
506	
507	# n = 2h + 2, m = 2h + 1: n & m is twice (h + 1) & h one width lower
508	def and_even(+w: Nat, +h: Nat) -> {SR.and_bits(1n+w, 2n+Nat.double(h), 1n+Nat.double(h)) == Nat.double(SR.and_bits(w, 1n+h, h)) : Nat}:
509	  %Equal.sym(Nat, C.bit(Nat.double(h)), 0n, bit_dbl0(h)) : {Nat.add(Nat.mul(_, C.bit(1n+Nat.double(h))), Nat.double(SR.and_bits(w, 1n+C.half(Nat.double(h)), C.half(1n+Nat.double(h))))) == Nat.double(SR.and_bits(w, 1n+h, h)) : Nat}
510	  %Equal.sym(Nat, C.half(Nat.double(h)), h, half_dbl0(h)) : {Nat.add(Nat.mul(0n, C.bit(1n+Nat.double(h))), Nat.double(SR.and_bits(w, 1n+_, C.half(1n+Nat.double(h))))) == Nat.double(SR.and_bits(w, 1n+h, h)) : Nat}
511	  %Equal.sym(Nat, C.half(1n+Nat.double(h)), h, half_dbl1(h)) : {Nat.add(Nat.mul(0n, C.bit(1n+Nat.double(h))), Nat.double(SR.and_bits(w, 1n+h, _))) == Nat.double(SR.and_bits(w, 1n+h, h)) : Nat}
512	  {==}
513	
514	# n = 2g + 3, m = 2g + 2: n & m is twice the low bits of g + 1
515	def and_odd(+w: Nat, +g: Nat) -> {SR.and_bits(1n+w, 3n+Nat.double(g), 2n+Nat.double(g)) == Nat.double(C.low(w, 1n+g)) : Nat}:
516	  %Equal.sym(Nat, C.bit(Nat.double(g)), 0n, bit_dbl0(g)) : {Nat.add(Nat.mul(C.bit(1n+Nat.double(g)), _), Nat.double(SR.and_bits(w, 1n+C.half(1n+Nat.double(g)), 1n+C.half(Nat.double(g))))) == Nat.double(C.low(w, 1n+g)) : Nat}
517	  %Equal.sym(Nat, C.half(Nat.double(g)), g, half_dbl0(g)) : {Nat.add(Nat.mul(C.bit(1n+Nat.double(g)), 0n), Nat.double(SR.and_bits(w, 1n+C.half(1n+Nat.double(g)), 1n+_))) == Nat.double(C.low(w, 1n+g)) : Nat}
518	  %Equal.sym(Nat, C.half(1n+Nat.double(g)), g, half_dbl1(g)) : {Nat.add(Nat.mul(C.bit(1n+Nat.double(g)), 0n), Nat.double(SR.and_bits(w, 1n+_, 1n+g))) == Nat.double(C.low(w, 1n+g)) : Nat}
519	  %Equal.sym(Nat, C.bit(1n+Nat.double(g)), 1n, bit_dbl1(g)) : {Nat.add(Nat.mul(_, 0n), Nat.double(SR.and_bits(w, 1n+g, 1n+g))) == Nat.double(C.low(w, 1n+g)) : Nat}
520	  Equal.cong(Nat, Nat, z => Nat.double(z), SR.and_bits(w, 1n+g, 1n+g), C.low(w, 1n+g), and_self(w, 1n+g))
521	
522	# THEOREM (the mask branch): with n = 1 + m, n & m == 0 and n < 2^v, each
523	# k < n is the low bits of floor(2^v / n) outputs
524	def mask(+v: Nat, +m: Nat, +k: Nat, +b: Nat, +h: Nat, +hb: {C.bit(m) == b : Nat}, +hh: {C.half(m) == h : Nat}, +hp: {SR.and_bits(v, 1n+m, m) == 0n : Nat}, +hw: {C.fits(v, 1n+m) == True{} : Bool}, +hk: {Nat.is_lt(k, 1n+m) == True{} : Bool}) -> {cm(v, m, k, C.shift(v, 1n)) == Nat.div(C.shift(v, 1n), 1n+m) : Nat}:
525	  match v:
526	    case 0n:
527	      Empty.absurd({cm(0n, m, k, 1n) == Nat.div(1n, 1n+m) : Nat}, L.false_true(hw))
528	    case 1n+w:
529	      match b h:
530	        case 0n 0n:
531	          +hm = decomp(m, 0n, 0n, hb, hh)
532	          %Equal.sym(Nat, m, 0n, hm) : {cm(1n+w, _, k, C.shift(1n+w, 1n)) == Nat.div(C.shift(1n+w, 1n), 1n+_) : Nat}
533	          Equal.trans(Nat, cm(1n+w, 0n, k, C.shift(1n+w, 1n)), C.shift(1n+w, 1n), Nat.div(C.shift(1n+w, 1n), 1n),
534	            cm0k(1n+w, k, C.shift(1n+w, 1n), Equal.trans(Bool, Nat.is_lt(k, 1n), Nat.is_lt(k, 1n+m), True{}, Equal.cong(Nat, Bool, z => Nat.is_lt(k, 1n+z), 0n, m, Equal.sym(Nat, m, 0n, hm)), hk)),
535	            Equal.sym(Nat, Nat.div(C.shift(1n+w, 1n), 1n), C.shift(1n+w, 1n), div_one(C.shift(1n+w, 1n))))
536	        case 0n 1n+g:
537	          +hm = decomp(m, 0n, 1n+g, hb, hh)
538	          +hp2 = Equal.trans(Nat, SR.and_bits(1n+w, 3n+Nat.double(g), 2n+Nat.double(g)), SR.and_bits(1n+w, 1n+m, m), 0n, Equal.cong(Nat, Nat, z => SR.and_bits(1n+w, 1n+z, z), 2n+Nat.double(g), m, Equal.sym(Nat, m, 2n+Nat.double(g), hm)), hp)
539	          +hw2 = Equal.trans(Bool, C.fits(w, 1n+g), C.fits(1n+w, 1n+m), True{},
540	            Equal.trans(Bool, C.fits(w, 1n+g), C.fits(w, 1n+C.half(1n+Nat.double(g))), C.fits(1n+w, 1n+m),
541	              Equal.cong(Nat, Bool, z => C.fits(w, 1n+z), g, C.half(1n+Nat.double(g)), Equal.sym(Nat, C.half(1n+Nat.double(g)), g, half_dbl1(g))),
542	              Equal.cong(Nat, Bool, z => C.fits(1n+w, 1n+z), 2n+Nat.double(g), m, Equal.sym(Nat, m, 2n+Nat.double(g), hm))),
543	            hw)
544	          +zero = Equal.trans(Nat, Nat.double(1n+g), Nat.double(C.low(w, 1n+g)), 0n,
545	            Equal.cong(Nat, Nat, z => Nat.double(z), 1n+g, C.low(w, 1n+g), Equal.sym(Nat, C.low(w, 1n+g), 1n+g, WI.low_fit(w, 1n+g, hw2))),
546	            Equal.trans(Nat, Nat.double(C.low(w, 1n+g)), SR.and_bits(1n+w, 3n+Nat.double(g), 2n+Nat.double(g)), 0n, Equal.sym(Nat, SR.and_bits(1n+w, 3n+Nat.double(g), 2n+Nat.double(g)), Nat.double(C.low(w, 1n+g)), and_odd(w, g)), hp2))
547	          Empty.absurd({cm(1n+w, m, k, C.shift(1n+w, 1n)) == Nat.div(C.shift(1n+w, 1n), 1n+m) : Nat}, N.succ_zero(g, dbl0(1n+g, zero)))
548	        case 1n _:
549	          +hm = decomp(m, 1n, h, hb, hh)
550	          +W = C.shift(w, 1n)
551	          +hp2 = Equal.trans(Nat, SR.and_bits(1n+w, 2n+Nat.double(h), 1n+Nat.double(h)), SR.and_bits(1n+w, 1n+m, m), 0n, Equal.cong(Nat, Nat, z => SR.and_bits(1n+w, 1n+z, z), 1n+Nat.double(h), m, Equal.sym(Nat, m, 1n+Nat.double(h), hm)), hp)
552	          +hp3 = dbl0(SR.and_bits(w, 1n+h, h), Equal.trans(Nat, Nat.double(SR.and_bits(w, 1n+h, h)), SR.and_bits(1n+w, 2n+Nat.double(h), 1n+Nat.double(h)), 0n, Equal.sym(Nat, SR.and_bits(1n+w, 2n+Nat.double(h), 1n+Nat.double(h)), Nat.double(SR.and_bits(w, 1n+h, h)), and_even(w, h)), hp2))
553	          +hw3 = Equal.trans(Bool, C.fits(w, 1n+h), C.fits(1n+w, 1n+m), True{},
554	            Equal.trans(Bool, C.fits(w, 1n+h), C.fits(w, 1n+C.half(Nat.double(h))), C.fits(1n+w, 1n+m),
555	              Equal.cong(Nat, Bool, z => C.fits(w, 1n+z), h, C.half(Nat.double(h)), Equal.sym(Nat, C.half(Nat.double(h)), h, half_dbl0(h))),
556	              Equal.cong(Nat, Bool, z => C.fits(1n+w, 1n+z), 1n+Nat.double(h), m, Equal.sym(Nat, m, 1n+Nat.double(h), hm))),
557	            hw)
558	          +hk3 = Equal.trans(Bool, Nat.is_lt(C.half(k), 1n+h), Nat.is_lt(k, 1n+m), True{},
559	            Equal.trans(Bool, Nat.is_lt(C.half(k), 1n+h), Nat.is_lt(k, Nat.double(1n+h)), Nat.is_lt(k, 1n+m), WI.lt_half(k, 1n+h),
560	              Equal.cong(Nat, Bool, z => Nat.is_lt(k, 1n+z), 1n+Nat.double(h), m, Equal.sym(Nat, m, 1n+Nat.double(h), hm))),
561	            hk)
562	          %Equal.sym(Nat, m, 1n+Nat.double(h), hm) : {cm(1n+w, _, k, Nat.double(W)) == Nat.div(Nat.double(W), 1n+_) : Nat}
563	          Equal.trans(Nat, cm(1n+w, 1n+Nat.double(h), k, Nat.double(W)), cm2(1n+w, 1n+Nat.double(h), k, W), Nat.div(Nat.double(W), 2n+Nat.double(h)),
564	            split(1n+w, 1n+Nat.double(h), k, W),
565	            Equal.trans(Nat, cm2(1n+w, 1n+Nat.double(h), k, W), cm(w, h, C.half(k), W), Nat.div(Nat.double(W), 2n+Nat.double(h)),
566	              pairs(w, h, k, W),
567	              Equal.trans(Nat, cm(w, h, C.half(k), W), Nat.div(W, 1n+h), Nat.div(Nat.double(W), 2n+Nat.double(h)),
568	                mask(w, h, C.half(k), C.bit(h), C.half(h), {==}, {==}, hp3, hw3, hk3),
569	                Equal.sym(Nat, Nat.div(Nat.double(W), 2n+Nat.double(h)), Nat.div(W, 1n+h), div_double(W, h)))))
570	        case 2n+q _:
571	          Empty.absurd({cm(1n+w, m, k, C.shift(1n+w, 1n)) == Nat.div(C.shift(1n+w, 1n), 1n+m) : Nat}, L.false_true(Equal.trans(Bool, Nat.is_le(2n+q, 1n), Nat.is_le(C.bit(m), 1n), True{}, Equal.cong(Nat, Bool, z => Nat.is_le(z, 1n), 2n+q, C.bit(m), Equal.sym(Nat, C.bit(m), 2n+q, hb)), WI.bit_le1(m))))
572	
573	# ---- both branches ----
574	
575	# the draw is the mask when n is a power of two
576	def draw_mask(+w: Nat, +m: Nat, +x: Nat, +hp: {Nat.is_eq(SR.and_bits(w, 1n+m, m), 0n) == True{} : Bool}) -> {SR.draw(w, x, 1n+m) == Some{SR.and_bits(w, x, m)} : Maybe<&2, Nat>}:
577	  %Equal.sym(Bool, Nat.is_eq(SR.and_bits(w, 1n+m, m), 0n), True{}, hp) : {SR.draw_pos(w, x, 1n+m, m, _) == Some{SR.and_bits(w, x, m)} : Maybe<&2, Nat>}
578	  {==}
579	
580	def count_mask(+w: Nat, +m: Nat, +k: Nat, +N: Nat, +hp: {Nat.is_eq(SR.and_bits(w, 1n+m, m), 0n) == True{} : Bool}) -> {SR.count(w, 1n+m, k, N) == cm(w, m, k, N) : Nat}:
581	  match N:
582	    case 0n:
583	      {==}
584	    case 1n+ +x:
585	      %Equal.sym(Nat, SR.count(w, 1n+m, k, x), cm(w, m, k, x), count_mask(w, m, k, x, hp)) : {Nat.add(_, SR.hit(SR.draw(w, x, 1n+m), k)) == Nat.add(cm(w, m, k, x), b2n(Nat.is_eq(SR.and_bits(w, x, m), k))) : Nat}
586	      %Equal.sym(Maybe<&2, Nat>, SR.draw(w, x, 1n+m), Some{SR.and_bits(w, x, m)}, draw_mask(w, m, x, hp)) : {Nat.add(cm(w, m, k, x), SR.hit(_, k)) == Nat.add(cm(w, m, k, x), b2n(Nat.is_eq(SR.and_bits(w, x, m), k))) : Nat}
587	      {==}
588	
589	def unbiased_c(+w: Nat, +m: Nat, +k: Nat, +hw: {C.fits(w, 1n+m) == True{} : Bool}, +hk: {Nat.is_lt(k, 1n+m) == True{} : Bool}, +c: Bool, +hc: {Nat.is_eq(SR.and_bits(w, 1n+m, m), 0n) == c : Bool}) -> {SR.count(w, 1n+m, k, C.shift(w, 1n)) == Nat.div(C.shift(w, 1n), 1n+m) : Nat}:
590	  match c:
591	    case False{}:
592	      lemire_count(w, m, k, hc, hw, hk)
593	    case True{}:
594	      Equal.trans(Nat, SR.count(w, 1n+m, k, C.shift(w, 1n)), cm(w, m, k, C.shift(w, 1n)), Nat.div(C.shift(w, 1n), 1n+m),
595	        count_mask(w, m, k, C.shift(w, 1n), hc),
596	        mask(w, m, k, C.bit(m), C.half(m), {==}, {==}, N.eq_from_is_eq(SR.and_bits(w, 1n+m, m), 0n, hc), hw, hk))
597	
598	# THEOREM (Lemire.unbiased): for every width w, bound 0 < n < 2^w and
599	# k < n, exactly floor(2^w / n) of the 2^w outputs draw k
600	def unbiased(+w: Nat, +n: Nat, +k: Nat, +hn: {Nat.is_lt(0n, n) == True{} : Bool}, +hw: {C.fits(w, n) == True{} : Bool}, +hk: {Nat.is_lt(k, n) == True{} : Bool}) -> {SR.count(w, n, k, C.shift(w, 1n)) == Nat.div(C.shift(w, 1n), n) : Nat}:
601	  match n:
602	    case 0n:
603	      Empty.absurd({SR.count(w, 0n, k, C.shift(w, 1n)) == Nat.div(C.shift(w, 1n), 0n) : Nat}, L.false_true(hn))
604	    case 1n+m:
605	      unbiased_c(w, m, k, hw, hk, Nat.is_eq(SR.and_bits(w, 1n+m, m), 0n), {==})