1	import Base
2	import ../../../spec/lib/common.bend as C
3	import ../../lib/nat.bend as N
4	import ../../lib/logic.bend as L
5	import ../../lib/lemmas/proofs/nat_algebra.bend as NA
6	import ../natural/sqrtn.bend as SQ
7	import ./width.bend as WW
8	import ./f64sqa.bend as AQ
9	import ./f64sqn.bend as QN
10	import ./w64dm.bend as DM
11	import ./w64div.bend as W64D
12	
13	# Nat facts for one step of Zimmermann's Karatsuba square root (used by
14	# f64sqr.bend): (a + b)^2 expanded, 2 (s 2^32) w == (w 2 s) 2^32,
15	# (s 2^32)^2 + r 2^64 == n 2^64, q <= 2^32, and (T + 3)^2 > (T + 1)^2 + 4 T.
16	
17	# (a + b)^2 expanded
18	def sqa(+a: Nat, +b: Nat) -> {Nat.mul(Nat.add(a, b), Nat.add(a, b)) == Nat.add(Nat.add(Nat.mul(a, a), Nat.add(Nat.mul(a, b), Nat.mul(a, b))), Nat.mul(b, b)) : Nat}:
19	  Equal.trans(Nat, Nat.mul(Nat.add(a, b), Nat.add(a, b)), Nat.add(Nat.mul(a, Nat.add(a, b)), Nat.mul(b, Nat.add(a, b))), Nat.add(Nat.add(Nat.mul(a, a), Nat.add(Nat.mul(a, b), Nat.mul(a, b))), Nat.mul(b, b)), NA.mul_add_right(a, b, Nat.add(a, b)), Equal.trans(Nat, Nat.add(Nat.mul(a, Nat.add(a, b)), Nat.mul(b, Nat.add(a, b))), Nat.add(Nat.add(Nat.mul(a, a), Nat.mul(a, b)), Nat.mul(b, Nat.add(a, b))), Nat.add(Nat.add(Nat.mul(a, a), Nat.add(Nat.mul(a, b), Nat.mul(a, b))), Nat.mul(b, b)), Equal.cong(Nat, Nat, z => Nat.add(z, Nat.mul(b, Nat.add(a, b))), Nat.mul(a, Nat.add(a, b)), Nat.add(Nat.mul(a, a), Nat.mul(a, b)), NA.mul_add_left(a, a, b)), Equal.trans(Nat, Nat.add(Nat.add(Nat.mul(a, a), Nat.mul(a, b)), Nat.mul(b, Nat.add(a, b))), Nat.add(Nat.add(Nat.mul(a, a), Nat.mul(a, b)), Nat.add(Nat.mul(b, a), Nat.mul(b, b))), Nat.add(Nat.add(Nat.mul(a, a), Nat.add(Nat.mul(a, b), Nat.mul(a, b))), Nat.mul(b, b)), Equal.cong(Nat, Nat, z => Nat.add(Nat.add(Nat.mul(a, a), Nat.mul(a, b)), z), Nat.mul(b, Nat.add(a, b)), Nat.add(Nat.mul(b, a), Nat.mul(b, b)), NA.mul_add_left(b, a, b)), Equal.trans(Nat, Nat.add(Nat.add(Nat.mul(a, a), Nat.mul(a, b)), Nat.add(Nat.mul(b, a), Nat.mul(b, b))), Nat.add(Nat.add(Nat.mul(a, a), Nat.mul(a, b)), Nat.add(Nat.mul(a, b), Nat.mul(b, b))), Nat.add(Nat.add(Nat.mul(a, a), Nat.add(Nat.mul(a, b), Nat.mul(a, b))), Nat.mul(b, b)), Equal.cong(Nat, Nat, z => Nat.add(Nat.add(Nat.mul(a, a), Nat.mul(a, b)), Nat.add(z, Nat.mul(b, b))), Nat.mul(b, a), Nat.mul(a, b), NA.mul_comm(b, a)), Equal.trans(Nat, Nat.add(Nat.add(Nat.mul(a, a), Nat.mul(a, b)), Nat.add(Nat.mul(a, b), Nat.mul(b, b))), Nat.add(Nat.mul(a, a), Nat.add(Nat.mul(a, b), Nat.add(Nat.mul(a, b), Nat.mul(b, b)))), Nat.add(Nat.add(Nat.mul(a, a), Nat.add(Nat.mul(a, b), Nat.mul(a, b))), Nat.mul(b, b)), NA.add_assoc(Nat.mul(a, a), Nat.mul(a, b), Nat.add(Nat.mul(a, b), Nat.mul(b, b))), Equal.trans(Nat, Nat.add(Nat.mul(a, a), Nat.add(Nat.mul(a, b), Nat.add(Nat.mul(a, b), Nat.mul(b, b)))), Nat.add(Nat.mul(a, a), 
20	
21	# 2 (s 2^32) w == (w D) 2^32 for D = 2 s
22	def uw2(+s: Nat, +w: Nat, +D: Nat, +hD: {Nat.add(s, s) == D : Nat}) -> {Nat.add(Nat.mul(C.shift(32n, s), w), Nat.mul(C.shift(32n, s), w)) == C.shift(32n, Nat.mul(w, D)) : Nat}:
23	  Equal.trans(Nat, Nat.add(Nat.mul(C.shift(32n, s), w), Nat.mul(C.shift(32n, s), w)), Nat.add(C.shift(32n, Nat.mul(s, w)), Nat.mul(C.shift(32n, s), w)), C.shift(32n, Nat.mul(w, D)), Equal.cong(Nat, Nat, z => Nat.add(z, Nat.mul(C.shift(32n, s), w)), Nat.mul(C.shift(32n, s), w), C.shift(32n, Nat.mul(s, w)), WW.shift_mul_l(32n, s, w)), Equal.trans(Nat, Nat.add(C.shift(32n, Nat.mul(s, w)), Nat.mul(C.shift(32n, s), w)), Nat.add(C.shift(32n, Nat.mul(s, w)), C.shift(32n, Nat.mul(s, w))), C.shift(32n, Nat.mul(w, D)), Equal.cong(Nat, Nat, z => Nat.add(C.shift(32n, Nat.mul(s, w)), z), Nat.mul(C.shift(32n, s), w), C.shift(32n, Nat.mul(s, w)), WW.shift_mul_l(32n, s, w)), Equal.trans(Nat, Nat.add(C.shift(32n, Nat.mul(s, w)), C.shift(32n, Nat.mul(s, w))), C.shift(32n, Nat.add(Nat.mul(s, w), Nat.mul(s, w))), C.shift(32n, Nat.mul(w, D)), Equal.sym(Nat, C.shift(32n, Nat.add(Nat.mul(s, w), Nat.mul(s, w))), Nat.add(C.shift(32n, Nat.mul(s, w)), C.shift(32n, Nat.mul(s, w))), WW.shift_add(32n, Nat.mul(s, w), Nat.mul(s, w))), Equal.trans(Nat, C.shift(32n, Nat.add(Nat.mul(s, w), Nat.mul(s, w))), C.shift(32n, Nat.mul(Nat.add(s, s), w)), C.shift(32n, Nat.mul(w, D)), Equal.cong(Nat, Nat, z => C.shift(32n, z), Nat.add(Nat.mul(s, w), Nat.mul(s, w)), Nat.mul(Nat.add(s, s), w), Equal.sym(Nat, Nat.mul(Nat.add(s, s), w), Nat.add(Nat.mul(s, w), Nat.mul(s, w)), NA.mul_add_right(s, s, w))), Equal.trans(Nat, C.shift(32n, Nat.mul(Nat.add(s, s), w)), C.shift(32n, Nat.mul(D, w)), C.shift(32n, Nat.mul(w, D)), Equal.cong(Nat, Nat, z => C.shift(32n, Nat.mul(z, w)), Nat.add(s, s), D, hD), Equal.cong(Nat, Nat, z => C.shift(32n, z), Nat.mul(D, w), Nat.mul(w, D), NA.mul_comm(D, w)))))))
24	
25	# (s 2^32)^2 + r 2^64 == n 2^64 for s^2 + r == n
26	def usq(+s: Nat, +r: Nat, +n: Nat, +hn: {Nat.add(Nat.mul(s, s), r) == n : Nat}) -> {Nat.add(Nat.mul(C.shift(32n, s), C.shift(32n, s)), C.shift(32n, C.shift(32n, r))) == C.shift(32n, C.shift(32n, n)) : Nat}:
27	  Equal.trans(Nat, Nat.add(Nat.mul(C.shift(32n, s), C.shift(32n, s)), C.shift(32n, C.shift(32n, r))), Nat.add(C.shift(32n, C.shift(32n, Nat.mul(s, s))), C.shift(32n, C.shift(32n, r))), C.shift(32n, C.shift(32n, n)), Equal.cong(Nat, Nat, z => Nat.add(z, C.shift(32n, C.shift(32n, r))), Nat.mul(C.shift(32n, s), C.shift(32n, s)), C.shift(32n, C.shift(32n, Nat.mul(s, s))), AQ.sh_sq(32n, s)), Equal.trans(Nat, Nat.add(C.shift(32n, C.shift(32n, Nat.mul(s, s))), C.shift(32n, C.shift(32n, r))), C.shift(32n, Nat.add(C.shift(32n, Nat.mul(s, s)), C.shift(32n, r))), C.shift(32n, C.shift(32n, n)), Equal.sym(Nat, C.shift(32n, Nat.add(C.shift(32n, Nat.mul(s, s)), C.shift(32n, r))), Nat.add(C.shift(32n, C.shift(32n, Nat.mul(s, s))), C.shift(32n, C.shift(32n, r))), WW.shift_add(32n, C.shift(32n, Nat.mul(s, s)), C.shift(32n, r))), Equal.trans(Nat, C.shift(32n, Nat.add(C.shift(32n, Nat.mul(s, s)), C.shift(32n, r))), C.shift(32n, C.shift(32n, Nat.add(Nat.mul(s, s), r))), C.shift(32n, C.shift(32n, n)), Equal.cong(Nat, Nat, z => C.shift(32n, z), Nat.add(C.shift(32n, Nat.mul(s, s)), C.shift(32n, r)), C.shift(32n, Nat.add(Nat.mul(s, s), r)), Equal.sym(Nat, C.shift(32n, Nat.add(Nat.mul(s, s), r)), Nat.add(C.shift(32n, Nat.mul(s, s)), C.shift(32n, r)), WW.shift_add(32n, Nat.mul(s, s), r))), Equal.cong(Nat, Nat, z => C.shift(32n, C.shift(32n, z)), Nat.add(Nat.mul(s, s), r), n, hn))))
28	
29	# from below: n 2^64 < (S0 + 1)^2, so isqrt(n 2^64) <= S0
30	# from above: S0^2 <= n 2^64 + q^2
31	# q <= 2^32: q D <= r 2^32 <= D 2^32 < (2^32 + 1) D
32	def qb(+one: Nat, +h1: {one == 1n : Nat}, +r: Nat, +q: Nat, +D: Nat, +hq: {Nat.is_le(Nat.mul(q, D), C.shift(32n, r)) == True{} : Bool}, +hr: {Nat.is_le(r, D) == True{} : Bool}, +hD: {Nat.is_lt(0n, D) == True{} : Bool}) -> {Nat.is_le(q, C.shift(32n, one)) == True{} : Bool}:
33	  +e1 = Equal.trans(Nat, C.shift(32n, D), Nat.mul(D, C.shift(32n, one)), Nat.mul(C.shift(32n, one), D), WW.shift_mul_one(32n, one, h1, D), NA.mul_comm(D, C.shift(32n, one)))
34	  +h2 = N.le_trans(Nat.mul(q, D), C.shift(32n, r), C.shift(32n, D), hq, WW.shift_mono(32n, r, D, hr))
35	  +h3 = L.subst(Nat, z => {Nat.is_le(Nat.mul(q, D), z) == True{} : Bool}, C.shift(32n, D), Nat.mul(C.shift(32n, one), D), e1, h2)
36	  +h4 = L.subst(Nat, z => {Nat.is_lt(Nat.mul(C.shift(32n, one), D), z) == True{} : Bool}, Nat.add(Nat.mul(C.shift(32n, one), D), D), Nat.add(D, Nat.mul(C.shift(32n, one), D)), NA.add_comm(Nat.mul(C.shift(32n, one), D), D), DM.lt_add_pos(Nat.mul(C.shift(32n, one), D), D, hD))
37	  +h5 = N.le_lt_trans(Nat.mul(q, D), Nat.mul(C.shift(32n, one), D), Nat.mul(1n+C.shift(32n, one), D), h3, h4)
38	  N.lt_succ_le(q, C.shift(32n, one), W64D.lt_cancel_mul(q, 1n+C.shift(32n, one), D, h5))
39	
40	def four(+x: Nat) -> {C.shift(2n, x) == Nat.add(Nat.add(x, x), Nat.add(x, x)) : Nat}:
41	  Equal.trans(Nat, Nat.double(Nat.double(x)), Nat.add(Nat.double(x), Nat.double(x)), Nat.add(Nat.add(x, x), Nat.add(x, x)), NA.double_self(Nat.double(x)), Equal.trans(Nat, Nat.add(Nat.double(x), Nat.double(x)), Nat.add(Nat.add(x, x), Nat.double(x)), Nat.add(Nat.add(x, x), Nat.add(x, x)), Equal.cong(Nat, Nat, z => Nat.add(z, Nat.double(x)), Nat.double(x), Nat.add(x, x), NA.double_self(x)), Equal.cong(Nat, Nat, z => Nat.add(Nat.add(x, x), z), Nat.double(x), Nat.add(x, x), NA.double_self(x))))
42	
43	def two_mul(+x: Nat) -> {Nat.mul(x, 2n) == Nat.add(x, x) : Nat}:
44	  QN.two_mul(x)
45	
46	# (T + 3)^2 exceeds (T + 1)^2 + 4 T
47	def t3(+T: Nat) -> {Nat.is_lt(Nat.add(Nat.mul(1n+T, 1n+T), Nat.add(Nat.add(T, T), Nat.add(T, T))), Nat.mul(Nat.add(1n+T, 2n), Nat.add(1n+T, 2n))) == True{} : Bool}:
48	  +a : Nat = 1n+T
49	  +e2 = Equal.trans(Nat, Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n)), Nat.add(Nat.add(a, a), Nat.mul(a, 2n)), Nat.add(Nat.add(a, a), Nat.add(a, a)), Equal.cong(Nat, Nat, z => Nat.add(z, Nat.mul(a, 2n)), Nat.mul(a, 2n), Nat.add(a, a), two_mul(a)), Equal.cong(Nat, Nat, z => Nat.add(Nat.add(a, a), z), Nat.mul(a, 2n), Nat.add(a, a), two_mul(a)))
50	  +hT = N.le_succ(T)
51	  +h4 = SQ.le_add2(Nat.add(T, T), Nat.add(a, a), Nat.add(T, T), Nat.add(a, a), SQ.le_add2(T, a, T, a, hT, hT), SQ.le_add2(T, a, T, a, hT, hT))
52	  +h5 = L.subst(Nat, z => {Nat.is_le(Nat.add(Nat.add(T, T), Nat.add(T, T)), z) == True{} : Bool}, Nat.add(Nat.add(a, a), Nat.add(a, a)), Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n)), Equal.sym(Nat, Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n)), Nat.add(Nat.add(a, a), Nat.add(a, a)), e2), h4)
53	  +h6 = N.le_lt_trans(Nat.add(Nat.add(T, T), Nat.add(T, T)), Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n)), Nat.add(Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n)), 4n), h5, DM.lt_add_pos(Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n)), 4n, {==}))
54	  +h7 = Equal.trans(Bool, Nat.is_lt(Nat.add(Nat.mul(a, a), Nat.add(Nat.add(T, T), Nat.add(T, T))), Nat.add(Nat.mul(a, a), Nat.add(Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n)), 4n))), Nat.is_lt(Nat.add(Nat.add(T, T), Nat.add(T, T)), Nat.add(Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n)), 4n)), True{}, WW.lt_cancel_l(Nat.mul(a, a), Nat.add(Nat.add(T, T), Nat.add(T, T)), Nat.add(Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n)), 4n)), h6)
55	  +h8 = L.subst(Nat, z => {Nat.is_lt(Nat.add(Nat.mul(a, a), Nat.add(Nat.add(T, T), Nat.add(T, T))), z) == True{} : Bool}, Nat.add(Nat.mul(a, a), Nat.add(Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n)), 4n)), Nat.add(Nat.add(Nat.mul(a, a), Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n))), Nat.mul(2n, 2n)), Equal.sym(Nat, Nat.add(Nat.add(Nat.mul(a, a), Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n))), Nat.mul(2n, 2n)), Nat.add(Nat.mul(a, a), Nat.add(Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n)), 4n)), NA.add_assoc(Nat.mul(a, a), Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n)), 4n)), h7)
56	  L.subst(Nat, z => {Nat.is_lt(Nat.add(Nat.mul(a, a), Nat.add(Nat.add(T, T), Nat.add(T, T))), z) == True{} : Bool}, Nat.add(Nat.add(Nat.mul(a, a), Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n))), Nat.mul(2n, 2n)), Nat.mul(Nat.add(a, 2n), Nat.add(a, 2n)), Equal.sym(Nat, Nat.mul(Nat.add(a, 2n), Nat.add(a, 2n)), Nat.add(Nat.add(Nat.mul(a, a), Nat.add(Nat.mul(a, 2n), Nat.mul(a, 2n))), Nat.mul(2n, 2n)), sqa(a, 2n)), h8)
57	
58	def k_hi_c(+X: Nat, +T: Nat, +S: Nat, +hlt: {Nat.is_lt(Nat.mul(S, S), Nat.add(Nat.mul(1n+T, 1n+T), Nat.add(Nat.add(T, T), Nat.add(T, T)))) == True{} : Bool}, +c: Bool, +hc: {Nat.is_le(S, Nat.add(T, 2n)) == c : Bool}) -> {c == True{} : Bool}:
59	  match c:
60	    case True{}:
61	      {==}
62	    case False{}:
63	      +h1 = N.lt_succ_le_succ(Nat.add(T, 2n), S, N.not_le_lt(S, Nat.add(T, 2n), hc))
64	      +h2 = AQ.sqmono(Nat.add(1n+T, 2n), S, h1)
65	      +h3 = N.lt_le_trans(Nat.mul(S, S), Nat.mul(Nat.add(1n+T, 2n), Nat.add(1n+T, 2n)), Nat.mul(S, S), N.lt_trans(Nat.mul(S, S), Nat.add(Nat.mul(1n+T, 1n+T), Nat.add(Nat.add(T, T), Nat.add(T, T))), Nat.mul(Nat.add(1n+T, 2n), Nat.add(1n+T, 2n)), hlt, t3(T)), h2)
66	      Empty.absurd({False{} == True{} : Bool}, DM.true_ne_false(Equal.trans(Bool, True{}, Nat.is_lt(Nat.mul(S, S), Nat.mul(S, S)), False{}, Equal.sym(Bool, Nat.is_lt(Nat.mul(S, S), Nat.mul(S, S)), True{}, h3), N.lt_irrefl(Nat.mul(S, S)))))
67	
68	# from above: S0 <= isqrt(n 2^64) + 2, for isqrt(n 2^64) >= 2^62