1	import Base
2	import ../../lib/common.bend as C
3	import ./field.bend as FS
4	
5	# The secp256k1 group (SEC 2 section 2.4.1): y^2 = x^3 + 7 over GF(p), the
6	# generator G and the order n, over the natural numbers (field.bend).
7	# Nothing of the implementation is used.
8	#
9	# Points are homogeneous projective triples (X : Y : Z) with x = X / Z,
10	# y = Y / Z and the point at infinity (0 : 1 : 0). Addition and doubling
11	# are Algorithms 7 and 9 of Renes, Costello and Batina, "Complete addition
12	# formulas for prime order elliptic curves" (EUROCRYPT 2016), for a = 0
13	# and b3 = 3 b = 21, written as the paper writes them: straight-line
14	# register programs (t0 <- X1 X2, ...), here with every assignment given a
15	# fresh register. This is HACL*'s Spec.K256.PointOps.point_add /
16	# point_double. Scalar multiplication is the textbook left-to-right
17	# double-and-add over the 256 bits of the scalar.
18	#
19	# That these projective formulas compute the textbook affine group law
20	# (SEC 1 section 2.2.1; Mathlib's WeierstrassCurve.Affine: slope, addX,
21	# addY) is not proved here: docs/CRYPTO_CONTRACTS.md lists it with the other
22	# group-law facts that are not proved.
23	
24	type SPoint is Data:
25	  SPoint{x: Nat, y: Nat, z: Nat}
26	
27	type SOp is Data:
28	  SAdd{i: Nat, j: Nat}
29	  SSub{i: Nat, j: Nat}
30	  SMul{i: Nat, j: Nat}
31	
32	def nth(i: Nat, regs: List<&2, Nat>) -> Nat:
33	  match i regs:
34	    case _ Nil{}: 0n
35	    case 0n r <> t: r
36	    case 1n+k r <> t: nth(k, t)
37	
38	def exec(+p: Nat, op: SOp, +regs: List<&2, Nat>) -> Nat:
39	  match op:
40	    case SAdd{i, j}: FS.madd(p, nth(i, regs), nth(j, regs))
41	    case SSub{i, j}: FS.msub(p, nth(i, regs), nth(j, regs))
42	    case SMul{i, j}: FS.mmul(p, nth(i, regs), nth(j, regs))
43	
44	# run a program: each instruction appends its result as a new register
45	def run_go(+p: Nat, prog: List<&2, SOp>, +regs: List<&2, Nat>) -> List<&2, Nat>:
46	  match prog:
47	    case Nil{}: regs
48	    case op <> t: run_go(p, t, List.append(&2, Nat, regs, [exec(p, op, regs)]))
49	
50	# (the first register is looked at first; the result is run_go's)
51	def run(+p: Nat, prog: List<&2, SOp>, regs: List<&2, Nat>) -> List<&2, Nat>:
52	  match regs:
53	    case Nil{}: run_go(p, prog, Nil{})
54	    case 0n <> t: run_go(p, prog, 0n <> t)
55	    case (1n+k) <> t: run_go(p, prog, (1n+k) <> t)
56	
57	# RCB Algorithm 7 (a = 0). Registers: 0 X1, 1 Y1, 2 Z1, 3 X2, 4 Y2, 5 Z2,
58	# 6 b3; the paper's lines, in order, write registers 7..39:
59	#   t0 = X1 X2 (7)      t1 = Y1 Y2 (8)      t2 = Z1 Z2 (9)
60	#   t3 = X1 + Y1 (10)   t4 = X2 + Y2 (11)   t3 = t3 t4 (12)
61	#   t4 = t0 + t1 (13)   t3 = t3 - t4 (14)   t4 = Y1 + Z1 (15)
62	#   X3 = Y2 + Z2 (16)   t4 = t4 X3 (17)     X3 = t1 + t2 (18)
63	#   t4 = t4 - X3 (19)   X3 = X1 + Z1 (20)   Y3 = X2 + Z2 (21)
64	#   X3 = X3 Y3 (22)     Y3 = t0 + t2 (23)   Y3 = X3 - Y3 (24)
65	#   X3 = t0 + t0 (25)   t0 = X3 + t0 (26)   t2 = b3 t2 (27)
66	#   Z3 = t1 + t2 (28)   t1 = t1 - t2 (29)   Y3 = b3 Y3 (30)
67	#   X3 = t4 Y3 (31)     t2 = t3 t1 (32)     X3 = t2 - X3 (33)
68	#   Y3 = Y3 t0 (34)     t1 = t1 Z3 (35)     Y3 = t1 + Y3 (36)
69	#   t0 = t0 t3 (37)     Z3 = Z3 t4 (38)     Z3 = Z3 + t0 (39)
70	# and the result is (X3 : Y3 : Z3) = registers (33 : 36 : 39).
71	def add_prog() -> List<&2, SOp>:
72	  [SMul{0n, 3n}, SMul{1n, 4n}, SMul{2n, 5n}, SAdd{0n, 1n}, SAdd{3n, 4n},
73	   SMul{10n, 11n}, SAdd{7n, 8n}, SSub{12n, 13n}, SAdd{1n, 2n}, SAdd{4n, 5n},
74	   SMul{15n, 16n}, SAdd{8n, 9n}, SSub{17n, 18n}, SAdd{0n, 2n}, SAdd{3n, 5n},
75	   SMul{20n, 21n}, SAdd{7n, 9n}, SSub{22n, 23n}, SAdd{7n, 7n}, SAdd{25n, 7n},
76	   SMul{6n, 9n}, SAdd{8n, 27n}, SSub{8n, 27n}, SMul{6n, 24n}, SMul{19n, 30n},
77	   SMul{14n, 29n}, SSub{32n, 31n}, SMul{30n, 26n}, SMul{29n, 28n}, SAdd{35n, 34n},
78	   SMul{26n, 14n}, SMul{28n, 19n}, SAdd{38n, 37n}]
79	
80	# RCB Algorithm 9 (a = 0). Registers: 0 X, 1 Y, 2 Z, 3 b3; lines write
81	# registers 4..21:
82	#   t0 = Y Y (4)        Z3 = t0 + t0 (5)    Z3 = Z3 + Z3 (6)
83	#   Z3 = Z3 + Z3 (7)    t1 = Y Z (8)        t2 = Z Z (9)
84	#   t2 = b3 t2 (10)     X3 = t2 Z3 (11)     Y3 = t0 + t2 (12)
85	#   Z3 = t1 Z3 (13)     t1 = t2 + t2 (14)   t2 = t1 + t2 (15)
86	#   t0 = t0 - t2 (16)   Y3 = t0 Y3 (17)     Y3 = X3 + Y3 (18)
87	#   t1 = X Y (19)       X3 = t0 t1 (20)     X3 = X3 + X3 (21)
88	# and the result is (X3 : Y3 : Z3) = registers (21 : 18 : 13).
89	def dbl_prog() -> List<&2, SOp>:
90	  [SMul{1n, 1n}, SAdd{4n, 4n}, SAdd{5n, 5n}, SAdd{6n, 6n}, SMul{1n, 2n},
91	   SMul{2n, 2n}, SMul{3n, 9n}, SMul{10n, 7n}, SAdd{4n, 10n}, SMul{8n, 7n},
92	   SAdd{10n, 10n}, SAdd{14n, 10n}, SSub{4n, 15n}, SMul{16n, 12n}, SAdd{11n, 17n},
93	   SMul{0n, 1n}, SMul{16n, 19n}, SAdd{20n, 20n}]
94	
95	def b3() -> Nat:
96	  21n
97	
98	def add_out(+r: List<&2, Nat>) -> SPoint:
99	  SPoint{nth(33n, r), nth(36n, r), nth(39n, r)}
100	
101	def padd(+p: Nat, a: SPoint, b: SPoint) -> SPoint:
102	  match a b:
103	    case SPoint{x1, y1, z1} SPoint{x2, y2, z2}:
104	      add_out(run(p, add_prog(), [x1, y1, z1, x2, y2, z2, b3()]))
105	
106	def dbl_out(+r: List<&2, Nat>) -> SPoint:
107	  SPoint{nth(21n, r), nth(18n, r), nth(13n, r)}
108	
109	def pdbl(+p: Nat, a: SPoint) -> SPoint:
110	  match a:
111	    case SPoint{x, y, z}:
112	      dbl_out(run(p, dbl_prog(), [x, y, z, b3()]))
113	
114	def infinity() -> SPoint:
115	  SPoint{0n, 1n, 0n}
116	
117	def is_inf(a: SPoint) -> Bool:
118	  match a:
119	    case SPoint{x, y, z}: Nat.is_eq(z, 0n)
120	
121	# bit i of k
122	def bit(+i: Nat, +k: Nat) -> Nat:
123	  C.bit(C.high(i, k))
124	
125	def step(b: Nat, +p: Nat, +a: SPoint, +d: SPoint) -> SPoint:
126	  match b:
127	    case 0n: d
128	    case 1n+c: padd(p, d, a)
129	
130	# bits i - 1, ..., 0 of k, most significant first (k is looked at first,
131	# so that on an unknown k nothing is unfolded)
132	def bits_go(i: Nat, +k: Nat) -> List<&2, Nat>:
133	  match i:
134	    case 0n: Nil{}
135	    case 1n+ +j: bit(j, k) <> bits_go(j, k)
136	
137	def bits(+i: Nat, k: Nat) -> List<&2, Nat>:
138	  match k:
139	    case 0n: bits_go(i, 0n)
140	    case 1n+j: bits_go(i, 1n+j)
141	
142	# left-to-right double-and-add: for each bit b, R = 2 R, then R = R + A
143	# when b is set
144	def lmul(+p: Nat, bs: List<&2, Nat>, +a: SPoint, r: SPoint) -> SPoint:
145	  match bs:
146	    case Nil{}: r
147	    case b <> t: lmul(p, t, a, step(b, p, a, pdbl(p, r)))
148	
149	# [k] A for 0 <= k < 2^256
150	def pmul(+p: Nat, +k: Nat, +a: SPoint) -> SPoint:
151	  lmul(p, bits(256n, k), a, infinity())
152	
153	# G (SEC 2), with z = 1
154	def g(+one: Nat) -> SPoint:
155	  SPoint{FS.digits(one, [6040n, 5880n, 33115n, 23026n, 10457n, 11726n, 64731n, 667n,
156	                         2823n, 52871n, 25237n, 21920n, 48044n, 63964n, 26238n, 31166n]),
157	         FS.digits(one, [54456n, 64272n, 53391n, 40007n, 21529n, 42629n, 46152n, 64791n,
158	                         2216n, 3601n, 64508n, 23972n, 50277n, 9891n, 55927n, 18490n]),
159	         one}
160	
161	# ---- affine coordinates ----
162	
163	type SAffine is Data:
164	  SAffine{x: Nat, y: Nat}
165	
166	def aff_z(+p: Nat, +x: Nat, +y: Nat, +zi: Nat) -> SAffine:
167	  SAffine{FS.mmul(p, x, zi), FS.mmul(p, y, zi)}
168	
169	# (X / Z, Y / Z), (0, 0) for the point at infinity
170	def to_affine(+p: Nat, a: SPoint) -> SAffine:
171	  match a:
172	    case SPoint{x, y, z}: aff_z(p, x, y, FS.minv(p, z))
173	
174	def aff_x(a: SAffine) -> Nat:
175	  match a:
176	    case SAffine{x, y}: x
177	
178	def aff_y(a: SAffine) -> Nat:
179	  match a:
180	    case SAffine{x, y}: y
181	
182	# x^3 + 7
183	def rhs(+p: Nat, +x: Nat) -> Nat:
184	  FS.madd(p, FS.mmul(p, FS.mmul(p, x, x), x), 7n)
185	
186	def on_curve(+p: Nat, +x: Nat, +y: Nat) -> Bool:
187	  Nat.is_eq(FS.mmul(p, y, y), rhs(p, x))
188	
189	# ---- SEC 1 section 2.3: octet strings ----
190	
191	def byte(+b: U32) -> Nat:
192	  U32.to_nat(U32.and(b, 255))
193	
194	# OS2IP (big-endian)
195	def os2ip_go(bs: List<&2, U32>, acc: Nat) -> Nat:
196	  match bs:
197	    case Nil{}: acc
198	    case b <> t: os2ip_go(t, Nat.add(C.shift(8n, acc), byte(b)))
199	
200	def os2ip(bs: List<&2, U32>) -> Nat:
201	  os2ip_go(bs, 0n)
202	
203	# I2OSP: the n low bytes of x, big-endian
204	def i2osp_le(n: Nat, +x: Nat) -> List<&2, U32>:
205	  match n:
206	    case 0n: Nil{}
207	    case 1n+k: U32.from_nat(C.low(8n, x)) <> i2osp_le(k, C.high(8n, x))
208	
209	def i2osp_go(n: Nat, +x: Nat) -> List<&2, U32>:
210	  List.reverse(&2, U32, i2osp_le(n, x))
211	
212	# (x is looked at first, so that on an unknown x nothing is unfolded)
213	def i2osp(+n: Nat, x: Nat) -> List<&2, U32>:
214	  match x:
215	    case 0n: i2osp_go(n, 0n)
216	    case 1n+k: i2osp_go(n, 1n+k)
217	
218	def length_is(n: Nat, bs: List<&2, U32>) -> Bool:
219	  Nat.is_eq(List.length(&2, U32, bs), n)
220	
221	def first(n: Nat, bs: List<&2, U32>) -> List<&2, U32>:
222	  match n bs:
223	    case 0n _: Nil{}
224	    case 1n+k Nil{}: Nil{}
225	    case 1n+k b <> t: b <> first(k, t)
226	
227	def after(n: Nat, bs: List<&2, U32>) -> List<&2, U32>:
228	  match n bs:
229	    case 0n _: bs
230	    case 1n+k Nil{}: Nil{}
231	    case 1n+k b <> t: after(k, t)
232	
233	def head(bs: List<&2, U32>) -> U32:
234	  match bs:
235	    case Nil{}: 0
236	    case b <> t: b
237	
238	# ---- SEC 1 section 2.3.3 / 2.3.4: point encodings ----
239	
240	def enc_c(a: SAffine) -> List<&2, U32>:
241	  match a:
242	    case SAffine{x, +y}: U32.from_nat(Nat.add(Nat.mod(y, 2n), 2n)) <> i2osp(32n, x)
243	
244	# compressed: 02 or 03 (the parity of y), then x
245	def encode_compressed(+p: Nat, a: SPoint) -> List<&2, U32>:
246	  enc_c(to_affine(p, a))
247	
248	def enc_u(a: SAffine) -> List<&2, U32>:
249	  match a:
250	    case SAffine{x, y}: 4 <> List.append(&2, U32, i2osp(32n, x), i2osp(32n, y))
251	
252	# uncompressed: 04, x, y
253	def encode_uncompressed(+p: Nat, a: SPoint) -> List<&2, U32>:
254	  enc_u(to_affine(p, a))
255	
256	# y = (x^3 + 7)^((p + 1) / 4), then y or p - y to get the wanted parity;
257	# no point when that y does not square to x^3 + 7
258	def pick_if(+p: Nat, +y: Nat, same: Bool) -> Nat:
259	  match same:
260	    case True{}: y
261	    case False{}: FS.mneg(p, y)
262	
263	def pick_parity(+p: Nat, +y: Nat, +par: Nat) -> Nat:
264	  pick_if(p, y, Nat.is_eq(Nat.mod(y, 2n), par))
265	
266	def decompress_if(+x: Nat, +y: Nat, ok: Bool) -> Maybe<&2, SPoint>:
267	  match ok:
268	    case True{}: Some{SPoint{x, y, 1n}}
269	    case False{}: None{}
270	
271	def decompress_y(+p: Nat, +x: Nat, +par: Nat, +y: Nat) -> Maybe<&2, SPoint>:
272	  decompress_if(x, pick_parity(p, y, par), Nat.is_eq(FS.mmul(p, y, y), rhs(p, x)))
273	
274	def decompress(+p: Nat, +x: Nat, +par: Nat) -> Maybe<&2, SPoint>:
275	  decompress_y(p, x, par, FS.fsqrt(p, rhs(p, x)))
276	
277	def decode_c_if(+p: Nat, +pre: U32, +x: Nat, ok: Bool) -> Maybe<&2, SPoint>:
278	  match ok:
279	    case True{}: decompress(p, x, U32.to_nat(U32.and(pre, 1)))
280	    case False{}: None{}
281	
282	def decode_c(+p: Nat, +pre: U32, +x: Nat) -> Maybe<&2, SPoint>:
283	  decode_c_if(p, pre, x, Bool.and(Bool.or(U32.is_eq(pre, 2), U32.is_eq(pre, 3)), Nat.is_lt(x, p)))
284	
285	def decode_u(+p: Nat, +pre: U32, +x: Nat, +y: Nat) -> Maybe<&2, SPoint>:
286	  decompress_if(x, y, Bool.and(Bool.and(U32.is_eq(pre, 4), Bool.and(Nat.is_lt(x, p), Nat.is_lt(y, p))), on_curve(p, x, y)))
287	
288	# a compressed (33 bytes) or uncompressed (65 bytes) public key; None when
289	# malformed, a coordinate is not below p, or the point is not on the curve
290	def decode_len(+p: Nat, +bs: List<&2, U32>, c33: Bool, c65: Bool) -> Maybe<&2, SPoint>:
291	  match c33 c65:
292	    case True{} _: decode_c(p, head(bs), os2ip(after(1n, bs)))
293	    case False{} True{}: decode_u(p, head(bs), os2ip(first(32n, after(1n, bs))), os2ip(after(33n, bs)))
294	    case False{} False{}: None{}
295	
296	def decode(+p: Nat, +bs: List<&2, U32>) -> Maybe<&2, SPoint>:
297	  decode_len(p, bs, length_is(33n, bs), length_is(65n, bs))