zzstoatzz.io
this repo has no description
1const std = @import("std");
2const Secp256k1 = std.crypto.ecc.Secp256k1;
3const Fe = @import("field.zig").Fe;
4const scalar = Secp256k1.scalar;
5
6const endo = @import("endo.zig");
7const point = @import("point.zig");
8const affine_mod = @import("affine.zig");
9const jacobian_mod = @import("jacobian.zig");
10const JacobianPoint = jacobian_mod.JacobianPoint;
11const AffinePoint = jacobian_mod.AffinePoint;
12
13pub const VerifyError = std.crypto.errors.IdentityElementError ||
14 std.crypto.errors.NonCanonicalError ||
15 error{SignatureVerificationFailed};
16
17/// Precomputed base point table: gTable()[i][j] = j * 256^i * G in affine.
18/// 32 subtables × 256 entries = 8192 affine points.
19/// Enables u1*G computation via byte-at-a-time lookup with zero doublings.
20/// Lazily initialized at runtime (u128 field arithmetic is too heavy
21/// for the comptime interpreter at this scale).
22var g_table_storage: [32][256]AffinePoint = undefined;
23var g_table_ready: bool = false;
24
25fn gTable() *const [32][256]AffinePoint {
26 if (!g_table_ready) {
27 g_table_storage = affine_mod.buildByteTable(Secp256k1.basePoint);
28 g_table_ready = true;
29 }
30 return &g_table_storage;
31}
32
33/// Reduce a 32-byte big-endian value to a scalar (mod n).
34fn reduceToScalar(h: [32]u8) scalar.Scalar {
35 var xs = [_]u8{0} ** 48;
36 @memcpy(xs[xs.len - 32 ..], &h);
37 return scalar.Scalar.fromBytes48(xs, .big);
38}
39
40/// Verify an ECDSA signature using precomputed tables and Jacobian arithmetic.
41///
42/// Optimizations:
43/// 1. u1*G via precomputed byte table: ~32 mixed adds, zero doublings
44/// 2. u2*Q via projective-table GLV: no field inversion, 4-bit windowed
45/// 3. Jacobian comparison (no field inversion for final check)
46/// 4. All field arithmetic via 5×52-bit limbs (libsecp256k1 style) — fast on 64-bit hardware
47pub fn verify(sig_r: [32]u8, sig_s: [32]u8, msg_hash: [32]u8, public_key: Secp256k1) VerifyError!void {
48 // parse and validate r, s
49 const r_sc = scalar.Scalar.fromBytes(sig_r, .big) catch return error.SignatureVerificationFailed;
50 const s_sc = scalar.Scalar.fromBytes(sig_s, .big) catch return error.SignatureVerificationFailed;
51 if (r_sc.isZero() or s_sc.isZero()) return error.IdentityElement;
52
53 // scalar_u1 = z * s^-1, scalar_u2 = r * s^-1
54 const z = reduceToScalar(msg_hash);
55 const s_inv = scalarInvert(s_sc);
56 const scalar_u1 = z.mul(s_inv).toBytes(.big);
57 const scalar_u2 = r_sc.mul(s_inv).toBytes(.little);
58
59 // 1. u1*G via precomputed byte table — direct 32-byte lookup, no GLV
60 const r1 = basePointMul(scalar_u1);
61
62 // 2. u2*Q via projective-table GLV — no field inversion
63 var split_u2 = endo.splitScalar(scalar_u2, .little) catch return error.SignatureVerificationFailed;
64 const zero_s = scalar.Scalar.zero.toBytes(.little);
65 var neg_p = false;
66 var neg_p_phi = false;
67 if (split_u2.r1[16] != 0) {
68 split_u2.r1 = scalar.neg(split_u2.r1, .little) catch zero_s;
69 neg_p = true;
70 }
71 if (split_u2.r2[16] != 0) {
72 split_u2.r2 = scalar.neg(split_u2.r2, .little) catch zero_s;
73 neg_p_phi = true;
74 }
75 const pk_affine26 = AffinePoint.fromStdlib(public_key.affineCoordinates());
76 const r2 = publicKeyMulProjective(split_u2, neg_p, neg_p_phi, pk_affine26);
77
78 // 3. combine results in Jacobian, compare without inversion
79 const q = r1.add(r2);
80
81 // reject identity
82 if (q.z.isZero()) return error.SignatureVerificationFailed;
83
84 // Jacobian comparison: r * Z² == X (mod p)
85 const r_fe = Fe.fromBytes(sig_r, .big);
86 if (!q.jacobianCompare(r_fe)) {
87 return error.SignatureVerificationFailed;
88 }
89}
90
91/// Compute u1*G using the precomputed byte table.
92/// Direct 32-byte lookup: G_TABLE[i][scalar[i]], no GLV decomposition.
93/// Cost: ~32 mixed Jacobian-affine additions, zero doublings.
94fn basePointMul(scalar_u1: [32]u8) JacobianPoint {
95 const table = gTable();
96 var acc = JacobianPoint.identity;
97
98 // table[i] corresponds to byte position i (little-endian).
99 // scalar_u1 is big-endian, so byte 31 is least significant.
100 // table[i] ↔ scalar_u1[31-i]
101 for (0..32) |i| {
102 const byte = scalar_u1[31 - i];
103 if (byte != 0) {
104 acc = acc.addMixed(table[i][byte]);
105 }
106 }
107
108 return acc;
109}
110
111/// Compute u2*Q using 2-way Jacobian Shamir with projective tables.
112/// Tables stay in Jacobian form — no batchToAffine field inversion.
113/// Cost: 128 Jacobian doublings + ~44 Jacobian additions.
114fn publicKeyMulProjective(
115 split: endo.SplitScalar,
116 neg_p: bool,
117 neg_p_phi: bool,
118 pk_affine: AffinePoint,
119) JacobianPoint {
120 // Build 9-entry Jacobian tables: [identity, 1P, 2P, ..., 8P]
121 const pk_table = point.precompute(pk_affine, 8);
122 const pk_phi_table = point.phiTableJacobian(9, pk_table);
123
124 const e1 = point.slide(split.r1);
125 const e2 = point.slide(split.r2);
126
127 var q = JacobianPoint.identity;
128 var pos: usize = 32; // 128-bit half-scalars → 32 nybbles + carry
129 while (true) : (pos -= 1) {
130 q = addSlotJ(q, &pk_table, e1[pos], neg_p);
131 q = addSlotJ(q, &pk_phi_table, e2[pos], neg_p_phi);
132 if (pos == 0) break;
133 q = q.dbl().dbl().dbl().dbl();
134 }
135 return q;
136}
137
138/// Add/subtract a Jacobian table entry based on signed digit and negation flag.
139/// Variable-time: branches on digit value (safe for public verification).
140inline fn addSlotJ(q: JacobianPoint, table: *const [9]JacobianPoint, slot: i8, negate: bool) JacobianPoint {
141 var s = slot;
142 if (negate) s = -s;
143 if (s > 0) {
144 return q.add(table[@intCast(s)]);
145 } else if (s < 0) {
146 return q.add(table[@intCast(-s)].negY());
147 }
148 return q;
149}
150
151/// Fermat scalar inversion: s^(n-2) mod n via addition chain.
152/// 253 squarings + 40 multiplications. Variable-time (safe for verification
153/// where all inputs are public).
154///
155/// Addition chain generated by github.com/mmcloughlin/addchain v0.4.0,
156/// ported from gitlab.com/yawning/secp256k1-voi.
157fn scalarInvert(x: scalar.Scalar) scalar.Scalar {
158 // Step 1: t0 = x^0x2
159 var t0 = x.sq();
160
161 // Step 2: t1 = x^0x3
162 var t1 = x.mul(t0);
163
164 // Step 3: t2 = x^0x5
165 var t2 = t0.mul(t1);
166
167 // Step 4: t3 = x^0x7
168 var t3 = t0.mul(t2);
169
170 // Step 5: t4 = x^0x9
171 var t4 = t0.mul(t3);
172
173 // Step 6: t5 = x^0xb
174 var t5 = t0.mul(t4);
175
176 // Step 7: t0 = x^0xd
177 t0 = t0.mul(t5);
178
179 // Step 9: t6 = x^0x34
180 var t6 = pow2k(t0, 2);
181
182 // Step 10: t6 = x^0x3f
183 t6 = t5.mul(t6);
184
185 // Step 11: t7 = x^0x7e
186 var t7 = t6.sq();
187
188 // Step 12: t7 = x^0x7f
189 t7 = x.mul(t7);
190
191 // Step 13: t8 = x^0xfe
192 var t8 = t7.sq();
193
194 // Step 14: t8 = x^0xff
195 t8 = x.mul(t8);
196
197 // Step 17: t9 = x^0x7f8
198 var t9 = pow2k(t8, 3);
199
200 // Step 19: t10 = x^0x1fe0
201 var t10 = pow2k(t9, 2);
202
203 // Step 20: t11 = x^0x3fc0
204 var t11 = t10.sq();
205
206 // Step 21: t12 = x^0x7f80
207 var t12 = t11.sq();
208
209 // Step 28: t13 = x^0x3fc000
210 const t13 = pow2k(t12, 7);
211
212 // Step 29: t11 = x^0x3fffc0
213 t11 = t11.mul(t13);
214
215 // Step 38: t11 = x^0x7fff8000
216 t11 = pow2k(t11, 9);
217
218 // Step 39: t12 = x^0x7fffff80
219 t12 = t12.mul(t11);
220
221 // Step 45: t11 = x^0x1fffffe000
222 t11 = pow2k(t12, 6);
223
224 // Step 46: t10 = x^0x1fffffffe0
225 t10 = t10.mul(t11);
226
227 // Step 72: t10 = x^0x7fffffff80000000
228 t10 = pow2k(t10, 26);
229
230 // Step 73: t12 = x^0x7fffffffffffff80
231 t12 = t12.mul(t10);
232
233 // Step 77: t10 = x^0x7fffffffffffff800
234 t10 = pow2k(t12, 4);
235
236 // Step 78: t9 = x^0x7fffffffffffffff8
237 t9 = t9.mul(t10);
238
239 // Step 138: t9 = x^0x7fffffffffffffff8000000000000000
240 t9 = pow2k(t9, 60);
241
242 // Step 139: t12 = x^0x7fffffffffffffffffffffffffffff80
243 t12 = t12.mul(t9);
244
245 // Step 140: t7 = x^0x7fffffffffffffffffffffffffffffff
246 t7 = t7.mul(t12);
247
248 // Step 145: t7 = x^0xfffffffffffffffffffffffffffffffe0
249 t7 = pow2k(t7, 5);
250
251 // Step 146: t7 = x^0xfffffffffffffffffffffffffffffffeb
252 t7 = t5.mul(t7);
253
254 // Step 149: t7 = x^0x7fffffffffffffffffffffffffffffff58
255 t7 = pow2k(t7, 3);
256
257 // Step 150: t7 = x^0x7fffffffffffffffffffffffffffffff5d
258 t7 = t2.mul(t7);
259
260 // Step 154: t7 = x^0x7fffffffffffffffffffffffffffffff5d0
261 t7 = pow2k(t7, 4);
262
263 // Step 155: t7 = x^0x7fffffffffffffffffffffffffffffff5d5
264 t7 = t2.mul(t7);
265
266 // Step 159: t7 = x^0x7fffffffffffffffffffffffffffffff5d50
267 t7 = pow2k(t7, 4);
268
269 // Step 160: t7 = x^0x7fffffffffffffffffffffffffffffff5d57
270 t7 = t3.mul(t7);
271
272 // Step 165: t7 = x^0xfffffffffffffffffffffffffffffffebaae0
273 t7 = pow2k(t7, 5);
274
275 // Step 166: t7 = x^0xfffffffffffffffffffffffffffffffebaaed
276 t7 = t0.mul(t7);
277
278 // Step 168: t7 = x^0x3fffffffffffffffffffffffffffffffaeabb4
279 t7 = pow2k(t7, 2);
280
281 // Step 169: t7 = x^0x3fffffffffffffffffffffffffffffffaeabb7
282 t7 = t1.mul(t7);
283
284 // Step 174: t7 = x^0x7fffffffffffffffffffffffffffffff5d576e0
285 t7 = pow2k(t7, 5);
286
287 // Step 175: t7 = x^0x7fffffffffffffffffffffffffffffff5d576e7
288 t7 = t3.mul(t7);
289
290 // Step 181: t7 = x^0x1fffffffffffffffffffffffffffffffd755db9c0
291 t7 = pow2k(t7, 6);
292
293 // Step 182: t7 = x^0x1fffffffffffffffffffffffffffffffd755db9cd
294 t7 = t0.mul(t7);
295
296 // Step 187: t7 = x^0x3fffffffffffffffffffffffffffffffaeabb739a0
297 t7 = pow2k(t7, 5);
298
299 // Step 188: t7 = x^0x3fffffffffffffffffffffffffffffffaeabb739ab
300 t7 = t5.mul(t7);
301
302 // Step 192: t7 = x^0x3fffffffffffffffffffffffffffffffaeabb739ab0
303 t7 = pow2k(t7, 4);
304
305 // Step 193: t7 = x^0x3fffffffffffffffffffffffffffffffaeabb739abd
306 t7 = t0.mul(t7);
307
308 // Step 196: t7 = x^0x1fffffffffffffffffffffffffffffffd755db9cd5e8
309 t7 = pow2k(t7, 3);
310
311 // Step 197: t7 = x^0x1fffffffffffffffffffffffffffffffd755db9cd5e9
312 t7 = x.mul(t7);
313
314 // Step 203: t7 = x^0x7fffffffffffffffffffffffffffffff5d576e7357a40
315 t7 = pow2k(t7, 6);
316
317 // Step 204: t2 = x^0x7fffffffffffffffffffffffffffffff5d576e7357a45
318 t2 = t2.mul(t7);
319
320 // Step 214: t2 = x^0x1fffffffffffffffffffffffffffffffd755db9cd5e91400
321 t2 = pow2k(t2, 10);
322
323 // Step 215: t2 = x^0x1fffffffffffffffffffffffffffffffd755db9cd5e91407
324 t2 = t3.mul(t2);
325
326 // Step 219: t2 = x^0x1fffffffffffffffffffffffffffffffd755db9cd5e914070
327 t2 = pow2k(t2, 4);
328
329 // Step 220: t3 = x^0x1fffffffffffffffffffffffffffffffd755db9cd5e914077
330 t3 = t3.mul(t2);
331
332 // Step 229: t3 = x^0x3fffffffffffffffffffffffffffffffaeabb739abd2280ee00
333 t3 = pow2k(t3, 9);
334
335 // Step 230: t8 = x^0x3fffffffffffffffffffffffffffffffaeabb739abd2280eeff
336 t8 = t8.mul(t3);
337
338 // Step 235: t8 = x^0x7fffffffffffffffffffffffffffffff5d576e7357a4501ddfe0
339 t8 = pow2k(t8, 5);
340
341 // Step 236: t8 = x^0x7fffffffffffffffffffffffffffffff5d576e7357a4501ddfe9
342 t8 = t4.mul(t8);
343
344 // Step 242: t8 = x^0x1fffffffffffffffffffffffffffffffd755db9cd5e9140777fa40
345 t8 = pow2k(t8, 6);
346
347 // Step 243: t5 = x^0x1fffffffffffffffffffffffffffffffd755db9cd5e9140777fa4b
348 t5 = t5.mul(t8);
349
350 // Step 247: t5 = x^0x1fffffffffffffffffffffffffffffffd755db9cd5e9140777fa4b0
351 t5 = pow2k(t5, 4);
352
353 // Step 248: t5 = x^0x1fffffffffffffffffffffffffffffffd755db9cd5e9140777fa4bd
354 t5 = t0.mul(t5);
355
356 // Step 253: t5 = x^0x3fffffffffffffffffffffffffffffffaeabb739abd2280eeff497a0
357 t5 = pow2k(t5, 5);
358
359 // Step 254: t1 = x^0x3fffffffffffffffffffffffffffffffaeabb739abd2280eeff497a3
360 t1 = t1.mul(t5);
361
362 // Step 260: t1 = x^0xfffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8c0
363 t1 = pow2k(t1, 6);
364
365 // Step 261: t1 = x^0xfffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd
366 t1 = t0.mul(t1);
367
368 // Step 271: t1 = x^0x3fffffffffffffffffffffffffffffffaeabb739abd2280eeff497a33400
369 t1 = pow2k(t1, 10);
370
371 // Step 272: t0 = x^0x3fffffffffffffffffffffffffffffffaeabb739abd2280eeff497a3340d
372 t0 = t0.mul(t1);
373
374 // Step 276: t0 = x^0x3fffffffffffffffffffffffffffffffaeabb739abd2280eeff497a3340d0
375 t0 = pow2k(t0, 4);
376
377 // Step 277: t4 = x^0x3fffffffffffffffffffffffffffffffaeabb739abd2280eeff497a3340d9
378 t4 = t4.mul(t0);
379
380 // Step 283: t4 = x^0xfffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd03640
381 t4 = pow2k(t4, 6);
382
383 // Step 284: t14 = x^0xfffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd03641
384 var t14 = x.mul(t4);
385
386 // Step 292: t14 = x^0xfffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364100
387 t14 = pow2k(t14, 8);
388
389 // Step 293: result = x^0xfffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd036413f
390 // This is x^(n-2) where n is the secp256k1 group order
391 return t6.mul(t14);
392}
393
394/// k repeated squarings: s^(2^k)
395inline fn pow2k(s: scalar.Scalar, comptime k: comptime_int) scalar.Scalar {
396 var r = s;
397 inline for (0..k) |_| {
398 r = r.sq();
399 }
400 return r;
401}
402
403test "scalarInvert matches stdlib invert" {
404 const S = scalar.Scalar;
405
406 // test with scalar = 1
407 const one = S.one;
408 const inv_one = scalarInvert(one);
409 const inv_one_std = one.invert();
410 try std.testing.expectEqual(inv_one.toBytes(.big), inv_one_std.toBytes(.big));
411
412 // test with a known scalar from base point
413 const s_bytes = Secp256k1.basePoint.affineCoordinates().x.toBytes(.big);
414 const s_val = S.fromBytes(s_bytes, .big) catch unreachable;
415 const inv_fermat = scalarInvert(s_val);
416 const inv_stdlib = s_val.invert();
417 try std.testing.expectEqual(inv_fermat.toBytes(.big), inv_stdlib.toBytes(.big));
418
419 // verify: s * s^-1 == 1
420 const product = s_val.mul(inv_fermat);
421 try std.testing.expectEqual(product.toBytes(.big), one.toBytes(.big));
422
423 // test with n-1 (the largest valid scalar)
424 // n-1 in big-endian
425 const n_minus_1_bytes = [_]u8{
426 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
427 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFE,
428 0xBA, 0xAE, 0xDC, 0xE6, 0xAF, 0x48, 0xA0, 0x3B,
429 0xBF, 0xD2, 0x5E, 0x8C, 0xD0, 0x36, 0x41, 0x40,
430 };
431 const nm1 = S.fromBytes(n_minus_1_bytes, .big) catch unreachable;
432 const inv_nm1_fermat = scalarInvert(nm1);
433 const inv_nm1_stdlib = nm1.invert();
434 try std.testing.expectEqual(inv_nm1_fermat.toBytes(.big), inv_nm1_stdlib.toBytes(.big));
435
436 // test with several more deterministic values
437 inline for ([_]u256{ 2, 7, 42, 0xDEADBEEF, 0x123456789ABCDEF0 }) |v| {
438 var buf: [32]u8 = undefined;
439 std.mem.writeInt(u256, &buf, v, .big);
440 const sv = S.fromBytes(buf, .big) catch unreachable;
441 const inv_f = scalarInvert(sv);
442 const inv_s = sv.invert();
443 try std.testing.expectEqual(inv_f.toBytes(.big), inv_s.toBytes(.big));
444 // verify round-trip
445 try std.testing.expectEqual(sv.mul(inv_f).toBytes(.big), one.toBytes(.big));
446 }
447}
448
449test "gTable()[0][1] is the base point" {
450 const G = Secp256k1.basePoint;
451 const g_affine = G.affineCoordinates();
452 const table = gTable();
453
454 const table_g = table[0][1];
455 const table_g_x = table_g.x.toStdlib();
456 const table_g_y = table_g.y.toStdlib();
457
458 try std.testing.expect(table_g_x.equivalent(g_affine.x));
459 try std.testing.expect(table_g_y.equivalent(g_affine.y));
460}
461
462test "gTable()[0][2] is 2G" {
463 const G = Secp256k1.basePoint;
464 const g2 = G.dbl().affineCoordinates();
465 const table = gTable();
466
467 const table_2g = table[0][2];
468 const table_2g_x = table_2g.x.toStdlib();
469 const table_2g_y = table_2g.y.toStdlib();
470
471 try std.testing.expect(table_2g_x.equivalent(g2.x));
472 try std.testing.expect(table_2g_y.equivalent(g2.y));
473}
474
475test "comptime addMixed gives correct 2G" {
476 @setEvalBranchQuota(100_000_000);
477 const g_affine = comptime AffinePoint.fromStdlib(Secp256k1.basePoint.affineCoordinates());
478 const j_g = comptime JacobianPoint.fromAffine(g_affine);
479 const j_2g = comptime j_g.addMixed(g_affine);
480
481 // convert to affine via toProjective
482 const proj = j_2g.toProjective();
483 const affine = proj.affineCoordinates();
484
485 const expected = Secp256k1.basePoint.dbl().affineCoordinates();
486
487 try std.testing.expect(affine.x.equivalent(expected.x));
488 try std.testing.expect(affine.y.equivalent(expected.y));
489}
490
491test "projective x-coordinate: X == x_affine * Z" {
492 // verify the coordinate convention: x_affine = X / Z (projective, not Jacobian)
493 const s = comptime blk: {
494 var buf: [32]u8 = undefined;
495 std.mem.writeInt(u256, &buf, 12345, .little);
496 break :blk buf;
497 };
498 const P = try Secp256k1.basePoint.mul(s, .little);
499 const affine = P.affineCoordinates();
500 const x_via_z = affine.x.mul(P.z);
501 try std.testing.expect(P.x.equivalent(x_via_z));
502}