1	import Base
2	import ./limbs.bend as L
3	import ./field.bend as F
4	import ./bytes.bend as B
5	
6	# The secp256k1 group y^2 = x^3 + 7 over GF(p) (SEC 2 section 2.4.1), in
7	# homogeneous projective coordinates (X : Y : Z), x = X / Z, y = Y / Z, the
8	# point at infinity (0 : 1 : 0). Addition and doubling are the complete
9	# formulas of Renes, Costello and Batina, "Complete addition formulas for
10	# prime order elliptic curves" (EUROCRYPT 2016), Algorithms 7 and 9 for
11	# a = 0 with b3 = 3 b = 21: one formula for every pair of inputs (equal
12	# points, opposite points and infinity included), so there is no
13	# exceptional case and no branch (the same choice as HACL*'s
14	# Hacl.Impl.K256.PointAdd / Spec.K256.PointOps).
15	#
16	# The formulas are straight-line programs over field registers, run by a
17	# small interpreter: registers 0.. are the inputs, each instruction appends
18	# one register. spec/crypto/secp256k1/curve.bend states the same programs,
19	# transcribed from the paper, over the field specification.
20	#
21	# Scalar multiplication is double-and-add over the 256 bits of the scalar,
22	# most significant first, adding in every step and selecting the result
23	# with an arithmetic mask (no branch and no table index on the scalar).
24	
25	type Op is Data:
26	  OAdd{i: Nat, j: Nat}
27	  OSub{i: Nat, j: Nat}
28	  OMul{i: Nat, j: Nat}
29	
30	type Point is Data:
31	  Point{x: List<&2, Nat>, y: List<&2, Nat>, z: List<&2, Nat>}
32	
33	def nth(i: Nat, regs: List<&2, List<&2, Nat>>) -> List<&2, Nat>:
34	  match i regs:
35	    case _ Nil{}: F.zero()
36	    case 0n r <> t: r
37	    case 1n+k r <> t: nth(k, t)
38	
39	def exec(op: Op, +regs: List<&2, List<&2, Nat>>) -> List<&2, Nat>:
40	  match op:
41	    case OAdd{i, j}: F.add(nth(i, regs), nth(j, regs))
42	    case OSub{i, j}: F.sub(nth(i, regs), nth(j, regs))
43	    case OMul{i, j}: F.mul(nth(i, regs), nth(j, regs))
44	
45	def run_go(prog: List<&2, Op>, +regs: List<&2, List<&2, Nat>>) -> List<&2, List<&2, Nat>>:
46	  match prog:
47	    case Nil{}: regs
48	    case op <> t: run_go(t, List.append(&2, List<&2, Nat>, regs, [exec(op, regs)]))
49	
50	# The first register is looked at before the program runs, so that the
51	# proof checker stops at once on unknown inputs (the result is run_go's).
52	def run(prog: List<&2, Op>, regs: List<&2, List<&2, Nat>>) -> List<&2, List<&2, Nat>>:
53	  match regs:
54	    case Nil{}: run_go(prog, Nil{})
55	    case Nil{} <> t: run_go(prog, Nil{} <> t)
56	    case (y <> u) <> t: run_go(prog, (y <> u) <> t)
57	
58	# RCB 2016 Algorithm 7 (a = 0): registers X1 Y1 Z1 X2 Y2 Z2 b3 = 0..6;
59	# the result is X3 Y3 Z3 = registers 33, 36, 39.
60	def add_prog() -> List<&2, Op>:
61	  [OMul{0n, 3n}, OMul{1n, 4n}, OMul{2n, 5n}, OAdd{0n, 1n}, OAdd{3n, 4n},
62	   OMul{10n, 11n}, OAdd{7n, 8n}, OSub{12n, 13n}, OAdd{1n, 2n}, OAdd{4n, 5n},
63	   OMul{15n, 16n}, OAdd{8n, 9n}, OSub{17n, 18n}, OAdd{0n, 2n}, OAdd{3n, 5n},
64	   OMul{20n, 21n}, OAdd{7n, 9n}, OSub{22n, 23n}, OAdd{7n, 7n}, OAdd{25n, 7n},
65	   OMul{6n, 9n}, OAdd{8n, 27n}, OSub{8n, 27n}, OMul{6n, 24n}, OMul{19n, 30n},
66	   OMul{14n, 29n}, OSub{32n, 31n}, OMul{30n, 26n}, OMul{29n, 28n}, OAdd{35n, 34n},
67	   OMul{26n, 14n}, OMul{28n, 19n}, OAdd{38n, 37n}]
68	
69	# RCB 2016 Algorithm 9 (a = 0): registers X Y Z b3 = 0..3; the result is
70	# X3 Y3 Z3 = registers 21, 18, 13.
71	def dbl_prog() -> List<&2, Op>:
72	  [OMul{1n, 1n}, OAdd{4n, 4n}, OAdd{5n, 5n}, OAdd{6n, 6n}, OMul{1n, 2n},
73	   OMul{2n, 2n}, OMul{3n, 9n}, OMul{10n, 7n}, OAdd{4n, 10n}, OMul{8n, 7n},
74	   OAdd{10n, 10n}, OAdd{14n, 10n}, OSub{4n, 15n}, OMul{16n, 12n}, OAdd{11n, 17n},
75	   OMul{0n, 1n}, OMul{16n, 19n}, OAdd{20n, 20n}]
76	
77	def b3() -> List<&2, Nat>:
78	  F.small(21n)
79	
80	def add_out(+r: List<&2, List<&2, Nat>>) -> Point:
81	  Point{nth(33n, r), nth(36n, r), nth(39n, r)}
82	
83	def add(p: Point, q: Point) -> Point:
84	  match p q:
85	    case Point{x1, y1, z1} Point{x2, y2, z2}:
86	      add_out(run(add_prog(), [x1, y1, z1, x2, y2, z2, b3()]))
87	
88	def dbl_out(+r: List<&2, List<&2, Nat>>) -> Point:
89	  Point{nth(21n, r), nth(18n, r), nth(13n, r)}
90	
91	def dbl(p: Point) -> Point:
92	  match p:
93	    case Point{x, y, z}:
94	      dbl_out(run(dbl_prog(), [x, y, z, b3()]))
95	
96	def infinity() -> Point:
97	  Point{F.zero(), F.one(), F.zero()}
98	
99	def neg(p: Point) -> Point:
100	  match p:
101	    case Point{x, y, z}: Point{x, F.neg(y), z}
102	
103	# b ? p : q, coordinate-wise and branch-free, for b in {0, 1}
104	def select(+b: Nat, p: Point, q: Point) -> Point:
105	  match p q:
106	    case Point{x1, y1, z1} Point{x2, y2, z2}:
107	      Point{F.select(b, x1, x2), F.select(b, y1, y2), F.select(b, z1, z2)}
108	
109	def step(+b: Nat, +p: Point, +d: Point) -> Point:
110	  select(b, add(d, p), d)
111	
112	def mul_go(bits: List<&2, Nat>, +p: Point, +r: Point) -> Point:
113	  match bits:
114	    case Nil{}: r
115	    case b <> t: mul_go(t, p, step(b, p, dbl(r)))
116	
117	# [k] p for a scalar k given by its 16 limbs (any value below 2^256)
118	def mul(k: List<&2, Nat>, +p: Point) -> Point:
119	  mul_go(L.bits(k), p, infinity())
120	
121	def gx() -> List<&2, Nat>:
122	  [6040n, 5880n, 33115n, 23026n, 10457n, 11726n, 64731n, 667n,
123	   2823n, 52871n, 25237n, 21920n, 48044n, 63964n, 26238n, 31166n]
124	
125	def gy() -> List<&2, Nat>:
126	  [54456n, 64272n, 53391n, 40007n, 21529n, 42629n, 46152n, 64791n,
127	   2216n, 3601n, 64508n, 23972n, 50277n, 9891n, 55927n, 18490n]
128	
129	# the generator G of SEC 2
130	def g() -> Point:
131	  Point{gx(), gy(), F.one()}
132	
133	def is_inf(p: Point) -> Bool:
134	  match p:
135	    case Point{x, y, z}: F.is_zero(z)
136	
137	# ---- affine coordinates and SEC 1 encodings ----
138	
139	type Affine is Data:
140	  Affine{x: List<&2, Nat>, y: List<&2, Nat>}
141	
142	def to_affine_z(+x: List<&2, Nat>, +y: List<&2, Nat>, +zi: List<&2, Nat>) -> Affine:
143	  Affine{F.mul(x, zi), F.mul(y, zi)}
144	
145	# (X / Z, Y / Z); (0, 0) for the point at infinity
146	def to_affine(p: Point) -> Affine:
147	  match p:
148	    case Point{x, y, z}: to_affine_z(x, y, F.inv(z))
149	
150	def aff_x(a: Affine) -> List<&2, Nat>:
151	  match a:
152	    case Affine{x, y}: x
153	
154	def aff_y(a: Affine) -> List<&2, Nat>:
155	  match a:
156	    case Affine{x, y}: y
157	
158	# x^3 + 7
159	def rhs(+x: List<&2, Nat>) -> List<&2, Nat>:
160	  F.add(F.mul(F.sq(x), x), F.small(7n))
161	
162	def on_curve(+x: List<&2, Nat>, +y: List<&2, Nat>) -> Bool:
163	  F.eq(F.sq(y), rhs(x))
164	
165	# the square root of x^3 + 7 with the given parity, if x^3 + 7 is a square
166	def lift_y(+y: List<&2, Nat>, +par: Nat) -> List<&2, Nat>:
167	  F.select(Nat.sub(Nat.add(F.parity(y), par), Nat.mul(Nat.mul(F.parity(y), par), 2n)), F.neg(y), y)
168	
169	def decompress_if(+x: List<&2, Nat>, +y: List<&2, Nat>, ok: Bool) -> Maybe<&2, Point>:
170	  match ok:
171	    case True{}: Some{Point{x, y, F.one()}}
172	    case False{}: None{}
173	
174	def decompress(+x: List<&2, Nat>, +par: Nat) -> Maybe<&2, Point>:
175	  +y = F.sqrt(rhs(x))
176	  decompress_if(x, lift_y(y, par), F.eq(F.sq(y), rhs(x)))
177	
178	def decode_c(+x: List<&2, Nat>, +par: Nat, ok: Bool) -> Maybe<&2, Point>:
179	  match ok:
180	    case True{}: decompress(x, par)
181	    case False{}: None{}
182	
183	def decode_u(+x: List<&2, Nat>, +y: List<&2, Nat>, ok: Bool) -> Maybe<&2, Point>:
184	  match ok:
185	    case True{}: Some{Point{x, y, F.one()}}
186	    case False{}: None{}
187	
188	def decode_33(+pre: U32, +rest: List<&2, U32>) -> Maybe<&2, Point>:
189	  +x = B.of_be(rest)
190	  decode_c(x, U32.to_nat(U32.and(pre, 1)), Bool.and(Bool.or(U32.is_eq(pre, 2), U32.is_eq(pre, 3)), F.lt_p(x)))
191	
192	def decode_65(+pre: U32, +rest: List<&2, U32>) -> Maybe<&2, Point>:
193	  +x = B.of_be(B.prefix(32n, rest))
194	  +y = B.of_be(B.suffix(32n, rest))
195	  decode_u(x, y, Bool.and(Bool.and(U32.is_eq(pre, 4), Bool.and(F.lt_p(x), F.lt_p(y))), on_curve(x, y)))
196	
197	def decode_len(+bs: List<&2, U32>, c33: Bool, c65: Bool) -> Maybe<&2, Point>:
198	  match c33 c65:
199	    case True{} _: decode_33(B.head(bs), B.tail(bs))
200	    case False{} True{}: decode_65(B.head(bs), B.tail(bs))
201	    case False{} False{}: None{}
202	
203	# SEC 1 section 2.3.4: a compressed (33 bytes, prefix 02/03) or uncompressed
204	# (65 bytes, prefix 04) public key, None when malformed, not below p, or
205	# not on the curve (the point at infinity, encoded 00, is refused)
206	def decode(+bs: List<&2, U32>) -> Maybe<&2, Point>:
207	  decode_len(bs, B.has_len(33n, bs), B.has_len(65n, bs))
208	
209	# SEC 1 section 2.3.3
210	def enc_c(a: Affine) -> List<&2, U32>:
211	  match a:
212	    case Affine{x, +y}: U32.from_nat(Nat.add(F.parity(y), 2n)) <> B.to_be(x)
213	
214	def encode_compressed(p: Point) -> List<&2, U32>:
215	  enc_c(to_affine(p))
216	
217	def enc_u(a: Affine) -> List<&2, U32>:
218	  match a:
219	    case Affine{x, y}: 4 <> List.append(&2, U32, B.to_be(x), B.to_be(y))
220	
221	def encode_uncompressed(p: Point) -> List<&2, U32>:
222	  enc_u(to_affine(p))