1	import Base
2	import ../lib/common.bend as C
3	import ../../src/math/natural.bend as M
4	import ../../src/math/number.bend as NB
5	
6	# Specification of src/math/number.bend: integer extras on Nat
7	# (python_style_math_stdlib_design.pdf sections 3.2, 8.7 and 10).
8	#
9	#   function    clauses                 statements mirror
10	#   bit_count   BitCount.value          Python int.bit_count; Mathlib
11	#                                       Nat.bits / List.count true
12	#   egcd        Egcd.gcd, Egcd.bezout   Knuth 4.5.2 Algorithm X; Mathlib
13	#                                       Nat.xgcd and Nat.gcd_eq_gcd_ab
14	#                                       (a * gcdA + b * gcdB = gcd)
15	#   is_prime    IsPrime.value           Mathlib Nat.Prime (2 <= n and no
16	#                                       divisor in [2, n)); the trial
17	#                                       division bound is Mathlib's
18	#                                       Nat.minFac_sq_le_self
19	
20	# the 1 bits among the k lowest bits of n
21	def ones(k: Nat, +n: Nat) -> Nat:
22	  match k:
23	    case 0n:
24	      0n
25	    case 1n+p:
26	      Nat.add(C.bit(n), ones(p, C.half(n)))
27	
28	# every n is below 2^n, so ones(n, n) counts all of n's 1 bits
29	def BitCount.value(+n: Nat) -> Type:
30	  {NB.bit_count(n) == ones(n, n) : Nat}
31	
32	# the gcd of the extended algorithm is the proved reference's
33	def eg_gcd(r: NB.EGcd) -> Nat:
34	  match r:
35	    case NB.EG{g, x, y, neg}:
36	      g
37	
38	def Egcd.gcd(+a: Nat, +b: Nat) -> Type:
39	  {eg_gcd(NB.egcd(a, b)) == M.gcd(a, b) : Nat}
40	
41	# the signed Bezout identity a*x - b*y == g (neg false) or b*y - a*x == g
42	# (neg true), stated without subtraction: the positive side equals g plus
43	# the negative side
44	def eg_pos(+a: Nat, +b: Nat, r: NB.EGcd) -> Nat:
45	  match r:
46	    case NB.EG{g, +x, +y, neg}:
47	      match neg:
48	        case True{}:
49	          Nat.mul(b, y)
50	        case False{}:
51	          Nat.mul(a, x)
52	
53	def eg_neg(+a: Nat, +b: Nat, r: NB.EGcd) -> Nat:
54	  match r:
55	    case NB.EG{+g, +x, +y, neg}:
56	      match neg:
57	        case True{}:
58	          Nat.add(g, Nat.mul(a, x))
59	        case False{}:
60	          Nat.add(g, Nat.mul(b, y))
61	
62	def Egcd.bezout(+a: Nat, +b: Nat) -> Type:
63	  {eg_pos(a, b, NB.egcd(a, b)) == eg_neg(a, b, NB.egcd(a, b)) : Nat}
64	
65	# no d in [d, d + k) divides n
66	def nodiv(k: Nat, +n: Nat, +d: Nat) -> Bool:
67	  match k:
68	    case 0n:
69	      True{}
70	    case 1n+p:
71	      Bool.and(Bool.not(Nat.is_eq(Nat.mod(n, d), 0n)), nodiv(p, n, 1n+d))
72	
73	# n is prime: 2 <= n and no d in [2, n) divides n
74	def prime(+n: Nat) -> Bool:
75	  match n:
76	    case 0n:
77	      False{}
78	    case 1n:
79	      False{}
80	    case 2n+ +p:
81	      nodiv(p, 2n+p, 2n)
82	
83	def IsPrime.value(+n: Nat) -> Type:
84	  {NB.is_prime(n) == prime(n) : Bool}