this repo has no description
feat: 10×26-bit field arithmetic (Fe26)
Replace stdlib 4×64-bit Montgomery field with direct 10×26-bit
representation for secp256k1. All point arithmetic, batch affine
conversion, endomorphism, and verification now operate in Fe26.
~9% faster than v0.0.1 baseline with safe normalize-on-output strategy.
Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
+1116
-181
+66
-39
src/affine.zig
+66
-39
src/affine.zig
···
1
1
const std = @import("std");
2
2
const Secp256k1 = std.crypto.ecc.Secp256k1;
3
-
const Fe = Secp256k1.Fe;
4
-
const AffineCoordinates = @TypeOf(Secp256k1.basePoint.affineCoordinates());
3
+
const Fe26 = @import("field.zig").Fe26;
4
+
const jacobian_mod = @import("jacobian.zig");
5
+
const JacobianPoint = jacobian_mod.JacobianPoint;
6
+
const AffinePoint26 = jacobian_mod.AffinePoint26;
5
7
6
-
/// Batch-convert projective points to affine using Montgomery's trick.
8
+
/// Batch-convert Jacobian points to affine using Montgomery's trick.
7
9
/// Uses 1 field inversion + 3*(n-1) field multiplications.
8
-
/// Identity points (z=0) are mapped to AffineCoordinates.identityElement.
9
-
pub fn batchToAffine(comptime n: usize, points: [n]Secp256k1) [n]AffineCoordinates {
10
+
/// For Jacobian: x_affine = X * Z^{-2}, y_affine = Y * Z^{-3}.
11
+
pub fn batchToAffine(comptime n: usize, points: [n]JacobianPoint) [n]AffinePoint26 {
10
12
// Replace zero Z with one to keep products invertible
11
-
var zs: [n]Fe = undefined;
13
+
var zs: [n]Fe26 = undefined;
12
14
for (0..n) |i| {
13
-
zs[i] = if (points[i].z.isZero()) Fe.one else points[i].z;
15
+
zs[i] = if (points[i].z.isZero()) Fe26.one else points[i].z;
14
16
}
15
17
16
18
// Forward: accumulate Z products
17
-
var products: [n]Fe = undefined;
19
+
var products: [n]Fe26 = undefined;
18
20
products[0] = zs[0];
19
21
for (1..n) |i| {
20
22
products[i] = products[i - 1].mul(zs[i]);
···
24
26
var inv = products[n - 1].invert();
25
27
26
28
// Backward: recover individual Z inverses and convert to affine
27
-
var result: [n]AffineCoordinates = undefined;
29
+
// For Jacobian: x = X * zinv², y = Y * zinv³
30
+
var result: [n]AffinePoint26 = undefined;
28
31
var i: usize = n - 1;
29
32
while (true) {
30
33
const z_inv = if (i > 0) inv.mul(products[i - 1]) else inv;
31
34
if (i > 0) inv = inv.mul(zs[i]);
32
35
33
-
result[i] = if (points[i].z.isZero())
34
-
AffineCoordinates.identityElement
35
-
else
36
-
.{ .x = points[i].x.mul(z_inv), .y = points[i].y.mul(z_inv) };
36
+
if (points[i].z.isZero()) {
37
+
result[i] = AffinePoint26.identity;
38
+
} else {
39
+
const z_inv2 = z_inv.sq();
40
+
const z_inv3 = z_inv2.mul(z_inv);
41
+
result[i] = .{
42
+
.x = points[i].x.mul(z_inv2),
43
+
.y = points[i].y.mul(z_inv3),
44
+
};
45
+
}
37
46
38
47
if (i == 0) break;
39
48
i -= 1;
···
43
52
}
44
53
45
54
/// Build a byte-indexed precomputed table for scalar multiplication.
46
-
/// table[i][j] = j * 256^i * base, stored as affine coordinates.
55
+
/// table[i][j] = j * 256^i * base, stored as AffinePoint26.
47
56
/// table[i][0] is the identity element (unused in lookups).
48
-
pub fn buildByteTable(base: Secp256k1) [16][256]AffineCoordinates {
49
-
// Phase 1: compute all 16*256 points in projective
50
-
var flat: [16 * 256]Secp256k1 = undefined;
57
+
pub fn buildByteTable(base: Secp256k1) [16][256]AffinePoint26 {
58
+
@setEvalBranchQuota(100_000_000);
51
59
52
-
var cur_base = base;
60
+
// Phase 1: compute all 16*256 points in Jacobian Fe26
61
+
var flat: [16 * 256]JacobianPoint = undefined;
62
+
63
+
var cur_base_affine = AffinePoint26.fromStdlib(base.affineCoordinates());
53
64
for (0..16) |sub| {
54
-
flat[sub * 256] = Secp256k1.identityElement;
55
-
flat[sub * 256 + 1] = cur_base;
65
+
flat[sub * 256] = JacobianPoint.identity;
66
+
flat[sub * 256 + 1] = JacobianPoint.fromAffine(cur_base_affine);
56
67
for (2..256) |j| {
57
-
flat[sub * 256 + j] = flat[sub * 256 + j - 1].add(cur_base);
68
+
flat[sub * 256 + j] = flat[sub * 256 + j - 1].addMixed(cur_base_affine);
58
69
}
59
70
if (sub < 15) {
60
71
// next subtable base = 256 * current base (8 doublings)
61
-
var next = cur_base;
72
+
var next = JacobianPoint.fromAffine(cur_base_affine);
62
73
next = next.dbl().dbl().dbl().dbl().dbl().dbl().dbl().dbl();
63
-
cur_base = next;
74
+
// convert back to affine for next iteration
75
+
const next_batch = batchToAffine(1, .{next});
76
+
cur_base_affine = next_batch[0];
64
77
}
65
78
}
66
79
···
68
81
const affine_flat = batchToAffine(16 * 256, flat);
69
82
70
83
// Reshape to [16][256]
71
-
var result: [16][256]AffineCoordinates = undefined;
84
+
var result: [16][256]AffinePoint26 = undefined;
72
85
for (0..16) |sub| {
73
86
for (0..256) |j| {
74
87
result[sub][j] = affine_flat[sub * 256 + j];
···
79
92
80
93
test "batchToAffine matches individual conversion" {
81
94
const G = Secp256k1.basePoint;
82
-
const points = [_]Secp256k1{
83
-
Secp256k1.identityElement,
84
-
G,
85
-
G.dbl(),
86
-
G.dbl().add(G),
87
-
};
95
+
const g_affine = G.affineCoordinates();
96
+
const g26 = AffinePoint26.fromStdlib(g_affine);
97
+
const StdAffineCoordinates = @TypeOf(g_affine);
98
+
99
+
const j_id = JacobianPoint.identity;
100
+
const j_g = JacobianPoint.fromAffine(g26);
101
+
const j_2g = j_g.dbl();
102
+
const j_3g = j_2g.addMixed(g26);
88
103
104
+
const points = [_]JacobianPoint{ j_id, j_g, j_2g, j_3g };
89
105
const result = batchToAffine(4, points);
90
106
91
107
// Identity
92
-
try std.testing.expect(result[0].x.equivalent(AffineCoordinates.identityElement.x));
108
+
try std.testing.expect(result[0].x.isZero());
93
109
94
-
// G
95
-
const g_affine = G.affineCoordinates();
96
-
try std.testing.expect(result[1].x.equivalent(g_affine.x));
97
-
try std.testing.expect(result[1].y.equivalent(g_affine.y));
110
+
// G — compare via stdlib
111
+
const result_g = StdAffineCoordinates{
112
+
.x = result[1].x.toStdlib(),
113
+
.y = result[1].y.toStdlib(),
114
+
};
115
+
try std.testing.expect(result_g.x.equivalent(g_affine.x));
116
+
try std.testing.expect(result_g.y.equivalent(g_affine.y));
98
117
99
118
// 2G
100
119
const g2_affine = G.dbl().affineCoordinates();
101
-
try std.testing.expect(result[2].x.equivalent(g2_affine.x));
102
-
try std.testing.expect(result[2].y.equivalent(g2_affine.y));
120
+
const result_2g = StdAffineCoordinates{
121
+
.x = result[2].x.toStdlib(),
122
+
.y = result[2].y.toStdlib(),
123
+
};
124
+
try std.testing.expect(result_2g.x.equivalent(g2_affine.x));
125
+
try std.testing.expect(result_2g.y.equivalent(g2_affine.y));
103
126
104
127
// 3G
105
128
const g3_affine = G.dbl().add(G).affineCoordinates();
106
-
try std.testing.expect(result[3].x.equivalent(g3_affine.x));
107
-
try std.testing.expect(result[3].y.equivalent(g3_affine.y));
129
+
const result_3g = StdAffineCoordinates{
130
+
.x = result[3].x.toStdlib(),
131
+
.y = result[3].y.toStdlib(),
132
+
};
133
+
try std.testing.expect(result_3g.x.equivalent(g3_affine.x));
134
+
try std.testing.expect(result_3g.y.equivalent(g3_affine.y));
108
135
}
+25
-14
src/endo.zig
+25
-14
src/endo.zig
···
1
1
const std = @import("std");
2
2
const Secp256k1 = std.crypto.ecc.Secp256k1;
3
-
const Fe = Secp256k1.Fe;
3
+
const Fe26 = @import("field.zig").Fe26;
4
+
const jacobian_mod = @import("jacobian.zig");
5
+
const AffinePoint26 = jacobian_mod.AffinePoint26;
6
+
const JacobianPoint = jacobian_mod.JacobianPoint;
4
7
5
8
/// Cube root of unity in GF(p): beta^3 = 1 mod p.
6
9
/// phi(x, y) = (beta * x, y) is equivalent to scalar multiplication by lambda,
7
10
/// but costs only 1 field multiply instead of ~65 doublings.
8
-
pub const beta = Fe.fromInt(
11
+
pub const beta = Fe26.fromInt(
9
12
55594575648329892869085402983802832744385952214688224221778511981742606582254,
10
-
) catch unreachable;
11
-
12
-
/// Apply the secp256k1 endomorphism: phi(P) = (beta * x, y, z).
13
-
/// Works in Jacobian coordinates — affine x is X/Z^2, so beta*(X/Z^2) = (beta*X)/Z^2.
14
-
pub fn phi(p: Secp256k1) Secp256k1 {
15
-
return .{ .x = p.x.mul(beta), .y = p.y, .z = p.z };
16
-
}
13
+
);
17
14
18
15
/// Apply the endomorphism to an affine point: (beta * x, y).
19
-
pub fn phiAffine(p: AffineCoordinates) AffineCoordinates {
16
+
pub fn phiAffine(p: AffinePoint26) AffinePoint26 {
20
17
return .{ .x = p.x.mul(beta), .y = p.y };
21
18
}
22
19
23
-
const AffineCoordinates = @TypeOf(Secp256k1.basePoint.affineCoordinates());
20
+
/// Apply the endomorphism to a Jacobian point: (beta * X, Y, Z).
21
+
/// phi preserves the Z coordinate since it only scales the x-coordinate.
22
+
pub fn phiJacobian(p: JacobianPoint) JacobianPoint {
23
+
return .{ .x = p.x.mul(beta), .y = p.y, .z = p.z };
24
+
}
24
25
25
26
/// Re-export stdlib's GLV scalar decomposition.
26
27
/// Decomposes k into (r1, r2) such that k = r1 + r2*lambda (mod n),
···
30
31
31
32
test "phi(G) matches lambda * G" {
32
33
const G = Secp256k1.basePoint;
33
-
const phi_G = phi(G);
34
+
const g_affine = G.affineCoordinates();
35
+
const g26 = AffinePoint26.fromStdlib(g_affine);
36
+
37
+
const phi_G = phiAffine(g26);
34
38
35
39
// lambda * G computed via scalar multiplication
36
40
const lambda_s = comptime s: {
···
44
48
break :s buf;
45
49
};
46
50
const lambda_G = try G.mulPublic(lambda_s, .little);
51
+
const lambda_G_affine = lambda_G.affineCoordinates();
47
52
48
-
// Both should represent the same point
49
-
try std.testing.expect(Secp256k1.equivalent(phi_G, lambda_G));
53
+
// compare via bytes
54
+
const phi_x = phi_G.x.normalize().toBytes(.big);
55
+
const lambda_x = Fe26.fromStdlib(lambda_G_affine.x).normalize().toBytes(.big);
56
+
try std.testing.expectEqualSlices(u8, &lambda_x, &phi_x);
57
+
58
+
const phi_y = phi_G.y.normalize().toBytes(.big);
59
+
const lambda_y = Fe26.fromStdlib(lambda_G_affine.y).normalize().toBytes(.big);
60
+
try std.testing.expectEqualSlices(u8, &lambda_y, &phi_y);
50
61
}
+717
src/field.zig
+717
src/field.zig
···
1
+
const std = @import("std");
2
+
const Secp256k1 = std.crypto.ecc.Secp256k1;
3
+
const StdFe = Secp256k1.Fe;
4
+
5
+
// field representation constants
6
+
const field_base = 26;
7
+
const field_words = 10;
8
+
const field_base_mask: u64 = (1 << field_base) - 1;
9
+
const field_msb_bits = 256 - (field_base * (field_words - 1));
10
+
const field_msb_mask: u64 = (1 << field_msb_bits) - 1;
11
+
12
+
// secp256k1 prime in 10x26 representation
13
+
const field_prime = [10]u32{
14
+
0x03fffc2f, 0x03ffffbf, 0x03ffffff, 0x03ffffff, 0x03ffffff,
15
+
0x03ffffff, 0x03ffffff, 0x03ffffff, 0x03ffffff, 0x003fffff,
16
+
};
17
+
18
+
/// Optimized secp256k1 field element: 10 × 26-bit limbs in 32-bit words.
19
+
///
20
+
/// 6 bits of overflow per word (10 in MSW) enables carry-free add/sub
21
+
/// for magnitudes up to 32. Mul/sq do inline reduction via the special
22
+
/// form p = 2^256 - 2^32 - 977, producing magnitude-1 output.
23
+
pub const Fe26 = struct {
24
+
n: [10]u32,
25
+
26
+
pub const zero: Fe26 = .{ .n = .{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } };
27
+
pub const one: Fe26 = .{ .n = .{ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 } };
28
+
29
+
/// Pack a 32-byte big-endian value into 10×26-bit limbs.
30
+
pub fn fromBytes(b: [32]u8, comptime endian: std.builtin.Endian) Fe26 {
31
+
const s = if (endian == .big) b else blk: {
32
+
var rev: [32]u8 = undefined;
33
+
for (0..32) |i| rev[i] = b[31 - i];
34
+
break :blk rev;
35
+
};
36
+
return .{ .n = .{
37
+
@as(u32, s[31]) | @as(u32, s[30]) << 8 | @as(u32, s[29]) << 16 |
38
+
(@as(u32, s[28]) & 0x3) << 24,
39
+
@as(u32, s[28]) >> 2 | @as(u32, s[27]) << 6 | @as(u32, s[26]) << 14 |
40
+
(@as(u32, s[25]) & 0xf) << 22,
41
+
@as(u32, s[25]) >> 4 | @as(u32, s[24]) << 4 | @as(u32, s[23]) << 12 |
42
+
(@as(u32, s[22]) & 0x3f) << 20,
43
+
@as(u32, s[22]) >> 6 | @as(u32, s[21]) << 2 | @as(u32, s[20]) << 10 |
44
+
@as(u32, s[19]) << 18,
45
+
@as(u32, s[18]) | @as(u32, s[17]) << 8 | @as(u32, s[16]) << 16 |
46
+
(@as(u32, s[15]) & 0x3) << 24,
47
+
@as(u32, s[15]) >> 2 | @as(u32, s[14]) << 6 | @as(u32, s[13]) << 14 |
48
+
(@as(u32, s[12]) & 0xf) << 22,
49
+
@as(u32, s[12]) >> 4 | @as(u32, s[11]) << 4 | @as(u32, s[10]) << 12 |
50
+
(@as(u32, s[9]) & 0x3f) << 20,
51
+
@as(u32, s[9]) >> 6 | @as(u32, s[8]) << 2 | @as(u32, s[7]) << 10 |
52
+
@as(u32, s[6]) << 18,
53
+
@as(u32, s[5]) | @as(u32, s[4]) << 8 | @as(u32, s[3]) << 16 |
54
+
(@as(u32, s[2]) & 0x3) << 24,
55
+
@as(u32, s[2]) >> 2 | @as(u32, s[1]) << 6 | @as(u32, s[0]) << 14,
56
+
} };
57
+
}
58
+
59
+
/// Unpack 10×26-bit limbs to a 32-byte big-endian value.
60
+
/// The field value MUST be normalized first.
61
+
pub fn toBytes(self: Fe26, comptime endian: std.builtin.Endian) [32]u8 {
62
+
var b: [32]u8 = undefined;
63
+
const n = self.n;
64
+
b[31] = @truncate(n[0]);
65
+
b[30] = @truncate(n[0] >> 8);
66
+
b[29] = @truncate(n[0] >> 16);
67
+
b[28] = @truncate((n[0] >> 24) & 0x3 | (n[1] & 0x3f) << 2);
68
+
b[27] = @truncate(n[1] >> 6);
69
+
b[26] = @truncate(n[1] >> 14);
70
+
b[25] = @truncate((n[1] >> 22) & 0xf | (n[2] & 0xf) << 4);
71
+
b[24] = @truncate(n[2] >> 4);
72
+
b[23] = @truncate(n[2] >> 12);
73
+
b[22] = @truncate((n[2] >> 20) & 0x3f | (n[3] & 0x3) << 6);
74
+
b[21] = @truncate(n[3] >> 2);
75
+
b[20] = @truncate(n[3] >> 10);
76
+
b[19] = @truncate(n[3] >> 18);
77
+
b[18] = @truncate(n[4]);
78
+
b[17] = @truncate(n[4] >> 8);
79
+
b[16] = @truncate(n[4] >> 16);
80
+
b[15] = @truncate((n[4] >> 24) & 0x3 | (n[5] & 0x3f) << 2);
81
+
b[14] = @truncate(n[5] >> 6);
82
+
b[13] = @truncate(n[5] >> 14);
83
+
b[12] = @truncate((n[5] >> 22) & 0xf | (n[6] & 0xf) << 4);
84
+
b[11] = @truncate(n[6] >> 4);
85
+
b[10] = @truncate(n[6] >> 12);
86
+
b[9] = @truncate((n[6] >> 20) & 0x3f | (n[7] & 0x3) << 6);
87
+
b[8] = @truncate(n[7] >> 2);
88
+
b[7] = @truncate(n[7] >> 10);
89
+
b[6] = @truncate(n[7] >> 18);
90
+
b[5] = @truncate(n[8]);
91
+
b[4] = @truncate(n[8] >> 8);
92
+
b[3] = @truncate(n[8] >> 16);
93
+
b[2] = @truncate((n[8] >> 24) & 0x3 | (n[9] & 0x3f) << 2);
94
+
b[1] = @truncate(n[9] >> 6);
95
+
b[0] = @truncate(n[9] >> 14);
96
+
if (endian == .little) {
97
+
var rev: [32]u8 = undefined;
98
+
for (0..32) |i| rev[i] = b[31 - i];
99
+
return rev;
100
+
}
101
+
return b;
102
+
}
103
+
104
+
/// Carry propagation + mod-p reduction. Two-pass: first propagate carries,
105
+
/// then conditionally subtract p.
106
+
pub fn normalize(self: Fe26) Fe26 {
107
+
// pass 1: reduce overflow from word 9 and propagate carries
108
+
var t9: u64 = self.n[9];
109
+
const m = t9 >> field_msb_bits;
110
+
t9 &= field_msb_mask;
111
+
var t0: u64 = self.n[0] + m * 977;
112
+
var t1: u64 = (t0 >> field_base) + self.n[1] + (m << 6);
113
+
t0 &= field_base_mask;
114
+
var t2: u64 = (t1 >> field_base) + self.n[2];
115
+
t1 &= field_base_mask;
116
+
var t3: u64 = (t2 >> field_base) + self.n[3];
117
+
t2 &= field_base_mask;
118
+
var t4: u64 = (t3 >> field_base) + self.n[4];
119
+
t3 &= field_base_mask;
120
+
var t5: u64 = (t4 >> field_base) + self.n[5];
121
+
t4 &= field_base_mask;
122
+
var t6: u64 = (t5 >> field_base) + self.n[6];
123
+
t5 &= field_base_mask;
124
+
var t7: u64 = (t6 >> field_base) + self.n[7];
125
+
t6 &= field_base_mask;
126
+
var t8: u64 = (t7 >> field_base) + self.n[8];
127
+
t7 &= field_base_mask;
128
+
t9 = (t8 >> field_base) + t9;
129
+
t8 &= field_base_mask;
130
+
131
+
// pass 2: conditional subtract p
132
+
// m=1 if value >= p, 0 otherwise
133
+
var m2 = @as(u64, @intFromBool(t9 == field_msb_mask));
134
+
m2 &= @intFromBool((t8 & t7 & t6 & t5 & t4 & t3 & t2) == field_base_mask);
135
+
m2 &= @intFromBool(t1 + 64 + ((t0 + 977) >> field_base) > field_base_mask);
136
+
m2 |= t9 >> field_msb_bits;
137
+
t0 += m2 * 977;
138
+
t1 = (t0 >> field_base) + t1 + (m2 << 6);
139
+
t0 &= field_base_mask;
140
+
t2 = (t1 >> field_base) + t2;
141
+
t1 &= field_base_mask;
142
+
t3 = (t2 >> field_base) + t3;
143
+
t2 &= field_base_mask;
144
+
t4 = (t3 >> field_base) + t4;
145
+
t3 &= field_base_mask;
146
+
t5 = (t4 >> field_base) + t5;
147
+
t4 &= field_base_mask;
148
+
t6 = (t5 >> field_base) + t6;
149
+
t5 &= field_base_mask;
150
+
t7 = (t6 >> field_base) + t7;
151
+
t6 &= field_base_mask;
152
+
t8 = (t7 >> field_base) + t8;
153
+
t7 &= field_base_mask;
154
+
t9 = (t8 >> field_base) + t9;
155
+
t8 &= field_base_mask;
156
+
t9 &= field_msb_mask;
157
+
158
+
return .{ .n = .{
159
+
@truncate(t0), @truncate(t1), @truncate(t2), @truncate(t3), @truncate(t4),
160
+
@truncate(t5), @truncate(t6), @truncate(t7), @truncate(t8), @truncate(t9),
161
+
} };
162
+
}
163
+
164
+
/// Returns true if the field value equals zero (after normalizing).
165
+
pub fn isZero(self: Fe26) bool {
166
+
const f = self.normalize();
167
+
const bits = f.n[0] | f.n[1] | f.n[2] | f.n[3] | f.n[4] |
168
+
f.n[5] | f.n[6] | f.n[7] | f.n[8] | f.n[9];
169
+
return bits == 0;
170
+
}
171
+
172
+
/// Returns true if two field values are equal (after normalizing both).
173
+
pub fn equivalent(self: Fe26, other: Fe26) bool {
174
+
const a = self.normalize();
175
+
const b = other.normalize();
176
+
const bits = (a.n[0] ^ b.n[0]) | (a.n[1] ^ b.n[1]) | (a.n[2] ^ b.n[2]) |
177
+
(a.n[3] ^ b.n[3]) | (a.n[4] ^ b.n[4]) | (a.n[5] ^ b.n[5]) |
178
+
(a.n[6] ^ b.n[6]) | (a.n[7] ^ b.n[7]) | (a.n[8] ^ b.n[8]) |
179
+
(a.n[9] ^ b.n[9]);
180
+
return bits == 0;
181
+
}
182
+
183
+
/// Carry-free addition. Magnitude increases by other's magnitude.
184
+
pub fn add(self: Fe26, other: Fe26) Fe26 {
185
+
return .{ .n = .{
186
+
self.n[0] + other.n[0], self.n[1] + other.n[1],
187
+
self.n[2] + other.n[2], self.n[3] + other.n[3],
188
+
self.n[4] + other.n[4], self.n[5] + other.n[5],
189
+
self.n[6] + other.n[6], self.n[7] + other.n[7],
190
+
self.n[8] + other.n[8], self.n[9] + other.n[9],
191
+
} };
192
+
}
193
+
194
+
/// Negate: result = (magnitude+1)*p - self. Magnitude increases by 1.
195
+
pub fn neg(self: Fe26, magnitude: u32) Fe26 {
196
+
const m = magnitude + 1;
197
+
return .{ .n = .{
198
+
m * field_prime[0] - self.n[0],
199
+
m * field_prime[1] - self.n[1],
200
+
m * field_prime[2] - self.n[2],
201
+
m * field_prime[3] - self.n[3],
202
+
m * field_prime[4] - self.n[4],
203
+
m * field_prime[5] - self.n[5],
204
+
m * field_prime[6] - self.n[6],
205
+
m * field_prime[7] - self.n[7],
206
+
m * field_prime[8] - self.n[8],
207
+
m * field_prime[9] - self.n[9],
208
+
} };
209
+
}
210
+
211
+
/// sub = self + (-other). Caller must provide the magnitude of `other` so
212
+
/// that neg produces non-negative limbs without needing normalize.
213
+
pub fn subMag(self: Fe26, other: Fe26, other_mag: u32) Fe26 {
214
+
return self.add(other.neg(other_mag));
215
+
}
216
+
217
+
/// sub with auto-normalize — safe but slower. Use subMag in hot paths.
218
+
pub fn sub(self: Fe26, other: Fe26) Fe26 {
219
+
return self.add(other.normalize().neg(1));
220
+
}
221
+
222
+
/// Double via carry-free addition.
223
+
pub fn dbl(self: Fe26) Fe26 {
224
+
return self.add(self);
225
+
}
226
+
227
+
/// Multiply by small integer. Magnitude scales by val.
228
+
pub fn mulInt(self: Fe26, val: u8) Fe26 {
229
+
const v: u32 = val;
230
+
return .{ .n = .{
231
+
self.n[0] * v, self.n[1] * v, self.n[2] * v, self.n[3] * v, self.n[4] * v,
232
+
self.n[5] * v, self.n[6] * v, self.n[7] * v, self.n[8] * v, self.n[9] * v,
233
+
} };
234
+
}
235
+
236
+
/// 100-product schoolbook multiplication with inline secp256k1 reduction.
237
+
/// Output magnitude: 1.
238
+
pub fn mul(self: Fe26, other: Fe26) Fe26 {
239
+
const a = self.n;
240
+
const b = other.n;
241
+
242
+
// schoolbook: 10×10 = terms for 2^(26*0) through 2^(26*18), plus overflow
243
+
var m: u64 = @as(u64, a[0]) * b[0];
244
+
const t0 = m & field_base_mask;
245
+
246
+
m = (m >> field_base) +
247
+
@as(u64, a[0]) * b[1] + @as(u64, a[1]) * b[0];
248
+
const t1 = m & field_base_mask;
249
+
250
+
m = (m >> field_base) +
251
+
@as(u64, a[0]) * b[2] + @as(u64, a[1]) * b[1] + @as(u64, a[2]) * b[0];
252
+
const t2 = m & field_base_mask;
253
+
254
+
m = (m >> field_base) +
255
+
@as(u64, a[0]) * b[3] + @as(u64, a[1]) * b[2] +
256
+
@as(u64, a[2]) * b[1] + @as(u64, a[3]) * b[0];
257
+
const t3 = m & field_base_mask;
258
+
259
+
m = (m >> field_base) +
260
+
@as(u64, a[0]) * b[4] + @as(u64, a[1]) * b[3] +
261
+
@as(u64, a[2]) * b[2] + @as(u64, a[3]) * b[1] + @as(u64, a[4]) * b[0];
262
+
const t4 = m & field_base_mask;
263
+
264
+
m = (m >> field_base) +
265
+
@as(u64, a[0]) * b[5] + @as(u64, a[1]) * b[4] +
266
+
@as(u64, a[2]) * b[3] + @as(u64, a[3]) * b[2] +
267
+
@as(u64, a[4]) * b[1] + @as(u64, a[5]) * b[0];
268
+
const t5 = m & field_base_mask;
269
+
270
+
m = (m >> field_base) +
271
+
@as(u64, a[0]) * b[6] + @as(u64, a[1]) * b[5] +
272
+
@as(u64, a[2]) * b[4] + @as(u64, a[3]) * b[3] +
273
+
@as(u64, a[4]) * b[2] + @as(u64, a[5]) * b[1] + @as(u64, a[6]) * b[0];
274
+
const t6 = m & field_base_mask;
275
+
276
+
m = (m >> field_base) +
277
+
@as(u64, a[0]) * b[7] + @as(u64, a[1]) * b[6] +
278
+
@as(u64, a[2]) * b[5] + @as(u64, a[3]) * b[4] +
279
+
@as(u64, a[4]) * b[3] + @as(u64, a[5]) * b[2] +
280
+
@as(u64, a[6]) * b[1] + @as(u64, a[7]) * b[0];
281
+
const t7 = m & field_base_mask;
282
+
283
+
m = (m >> field_base) +
284
+
@as(u64, a[0]) * b[8] + @as(u64, a[1]) * b[7] +
285
+
@as(u64, a[2]) * b[6] + @as(u64, a[3]) * b[5] +
286
+
@as(u64, a[4]) * b[4] + @as(u64, a[5]) * b[3] +
287
+
@as(u64, a[6]) * b[2] + @as(u64, a[7]) * b[1] + @as(u64, a[8]) * b[0];
288
+
const t8 = m & field_base_mask;
289
+
290
+
m = (m >> field_base) +
291
+
@as(u64, a[0]) * b[9] + @as(u64, a[1]) * b[8] +
292
+
@as(u64, a[2]) * b[7] + @as(u64, a[3]) * b[6] +
293
+
@as(u64, a[4]) * b[5] + @as(u64, a[5]) * b[4] +
294
+
@as(u64, a[6]) * b[3] + @as(u64, a[7]) * b[2] +
295
+
@as(u64, a[8]) * b[1] + @as(u64, a[9]) * b[0];
296
+
const t9 = m & field_base_mask;
297
+
298
+
m = (m >> field_base) +
299
+
@as(u64, a[1]) * b[9] + @as(u64, a[2]) * b[8] +
300
+
@as(u64, a[3]) * b[7] + @as(u64, a[4]) * b[6] +
301
+
@as(u64, a[5]) * b[5] + @as(u64, a[6]) * b[4] +
302
+
@as(u64, a[7]) * b[3] + @as(u64, a[8]) * b[2] + @as(u64, a[9]) * b[1];
303
+
const t10 = m & field_base_mask;
304
+
305
+
m = (m >> field_base) +
306
+
@as(u64, a[2]) * b[9] + @as(u64, a[3]) * b[8] +
307
+
@as(u64, a[4]) * b[7] + @as(u64, a[5]) * b[6] +
308
+
@as(u64, a[6]) * b[5] + @as(u64, a[7]) * b[4] +
309
+
@as(u64, a[8]) * b[3] + @as(u64, a[9]) * b[2];
310
+
const t11 = m & field_base_mask;
311
+
312
+
m = (m >> field_base) +
313
+
@as(u64, a[3]) * b[9] + @as(u64, a[4]) * b[8] +
314
+
@as(u64, a[5]) * b[7] + @as(u64, a[6]) * b[6] +
315
+
@as(u64, a[7]) * b[5] + @as(u64, a[8]) * b[4] + @as(u64, a[9]) * b[3];
316
+
const t12 = m & field_base_mask;
317
+
318
+
m = (m >> field_base) +
319
+
@as(u64, a[4]) * b[9] + @as(u64, a[5]) * b[8] +
320
+
@as(u64, a[6]) * b[7] + @as(u64, a[7]) * b[6] +
321
+
@as(u64, a[8]) * b[5] + @as(u64, a[9]) * b[4];
322
+
const t13 = m & field_base_mask;
323
+
324
+
m = (m >> field_base) +
325
+
@as(u64, a[5]) * b[9] + @as(u64, a[6]) * b[8] +
326
+
@as(u64, a[7]) * b[7] + @as(u64, a[8]) * b[6] + @as(u64, a[9]) * b[5];
327
+
const t14 = m & field_base_mask;
328
+
329
+
m = (m >> field_base) +
330
+
@as(u64, a[6]) * b[9] + @as(u64, a[7]) * b[8] +
331
+
@as(u64, a[8]) * b[7] + @as(u64, a[9]) * b[6];
332
+
const t15 = m & field_base_mask;
333
+
334
+
m = (m >> field_base) +
335
+
@as(u64, a[7]) * b[9] + @as(u64, a[8]) * b[8] + @as(u64, a[9]) * b[7];
336
+
const t16 = m & field_base_mask;
337
+
338
+
m = (m >> field_base) +
339
+
@as(u64, a[8]) * b[9] + @as(u64, a[9]) * b[8];
340
+
const t17 = m & field_base_mask;
341
+
342
+
m = (m >> field_base) + @as(u64, a[9]) * b[9];
343
+
const t18 = m & field_base_mask;
344
+
345
+
const t19 = m >> field_base;
346
+
347
+
return reduce(t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17, t18, t19);
348
+
}
349
+
350
+
/// 55-unique-product symmetric square with inline secp256k1 reduction.
351
+
/// Exploits a[i]*a[j] == a[j]*a[i] symmetry. Output magnitude: 1.
352
+
pub fn sq(self: Fe26) Fe26 {
353
+
const a = self.n;
354
+
355
+
var m: u64 = @as(u64, a[0]) * a[0];
356
+
const t0 = m & field_base_mask;
357
+
358
+
m = (m >> field_base) + 2 * @as(u64, a[0]) * a[1];
359
+
const t1 = m & field_base_mask;
360
+
361
+
m = (m >> field_base) +
362
+
2 * @as(u64, a[0]) * a[2] + @as(u64, a[1]) * a[1];
363
+
const t2 = m & field_base_mask;
364
+
365
+
m = (m >> field_base) +
366
+
2 * @as(u64, a[0]) * a[3] + 2 * @as(u64, a[1]) * a[2];
367
+
const t3 = m & field_base_mask;
368
+
369
+
m = (m >> field_base) +
370
+
2 * @as(u64, a[0]) * a[4] + 2 * @as(u64, a[1]) * a[3] +
371
+
@as(u64, a[2]) * a[2];
372
+
const t4 = m & field_base_mask;
373
+
374
+
m = (m >> field_base) +
375
+
2 * @as(u64, a[0]) * a[5] + 2 * @as(u64, a[1]) * a[4] +
376
+
2 * @as(u64, a[2]) * a[3];
377
+
const t5 = m & field_base_mask;
378
+
379
+
m = (m >> field_base) +
380
+
2 * @as(u64, a[0]) * a[6] + 2 * @as(u64, a[1]) * a[5] +
381
+
2 * @as(u64, a[2]) * a[4] + @as(u64, a[3]) * a[3];
382
+
const t6 = m & field_base_mask;
383
+
384
+
m = (m >> field_base) +
385
+
2 * @as(u64, a[0]) * a[7] + 2 * @as(u64, a[1]) * a[6] +
386
+
2 * @as(u64, a[2]) * a[5] + 2 * @as(u64, a[3]) * a[4];
387
+
const t7 = m & field_base_mask;
388
+
389
+
m = (m >> field_base) +
390
+
2 * @as(u64, a[0]) * a[8] + 2 * @as(u64, a[1]) * a[7] +
391
+
2 * @as(u64, a[2]) * a[6] + 2 * @as(u64, a[3]) * a[5] +
392
+
@as(u64, a[4]) * a[4];
393
+
const t8 = m & field_base_mask;
394
+
395
+
m = (m >> field_base) +
396
+
2 * @as(u64, a[0]) * a[9] + 2 * @as(u64, a[1]) * a[8] +
397
+
2 * @as(u64, a[2]) * a[7] + 2 * @as(u64, a[3]) * a[6] +
398
+
2 * @as(u64, a[4]) * a[5];
399
+
const t9 = m & field_base_mask;
400
+
401
+
m = (m >> field_base) +
402
+
2 * @as(u64, a[1]) * a[9] + 2 * @as(u64, a[2]) * a[8] +
403
+
2 * @as(u64, a[3]) * a[7] + 2 * @as(u64, a[4]) * a[6] +
404
+
@as(u64, a[5]) * a[5];
405
+
const t10 = m & field_base_mask;
406
+
407
+
m = (m >> field_base) +
408
+
2 * @as(u64, a[2]) * a[9] + 2 * @as(u64, a[3]) * a[8] +
409
+
2 * @as(u64, a[4]) * a[7] + 2 * @as(u64, a[5]) * a[6];
410
+
const t11 = m & field_base_mask;
411
+
412
+
m = (m >> field_base) +
413
+
2 * @as(u64, a[3]) * a[9] + 2 * @as(u64, a[4]) * a[8] +
414
+
2 * @as(u64, a[5]) * a[7] + @as(u64, a[6]) * a[6];
415
+
const t12 = m & field_base_mask;
416
+
417
+
m = (m >> field_base) +
418
+
2 * @as(u64, a[4]) * a[9] + 2 * @as(u64, a[5]) * a[8] +
419
+
2 * @as(u64, a[6]) * a[7];
420
+
const t13 = m & field_base_mask;
421
+
422
+
m = (m >> field_base) +
423
+
2 * @as(u64, a[5]) * a[9] + 2 * @as(u64, a[6]) * a[8] +
424
+
@as(u64, a[7]) * a[7];
425
+
const t14 = m & field_base_mask;
426
+
427
+
m = (m >> field_base) +
428
+
2 * @as(u64, a[6]) * a[9] + 2 * @as(u64, a[7]) * a[8];
429
+
const t15 = m & field_base_mask;
430
+
431
+
m = (m >> field_base) +
432
+
2 * @as(u64, a[7]) * a[9] + @as(u64, a[8]) * a[8];
433
+
const t16 = m & field_base_mask;
434
+
435
+
m = (m >> field_base) + 2 * @as(u64, a[8]) * a[9];
436
+
const t17 = m & field_base_mask;
437
+
438
+
m = (m >> field_base) + @as(u64, a[9]) * a[9];
439
+
const t18 = m & field_base_mask;
440
+
441
+
const t19 = m >> field_base;
442
+
443
+
return reduce(t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17, t18, t19);
444
+
}
445
+
446
+
/// Inline secp256k1 reduction of a 20-word product.
447
+
/// p = 2^256 - c where c = 4294968273 = 2^26 * 64 + 977.
448
+
/// Upper terms (t10..t19) at bit 260+, so adjusted c = c * 16:
449
+
/// c_lo = 977 * 16 = 15632, c_hi = 64 * 16 = 1024.
450
+
/// Final term t19 uses full c * 16^2 = 68719492368.
451
+
fn reduce(
452
+
t0_in: u64,
453
+
t1_in: u64,
454
+
t2_in: u64,
455
+
t3_in: u64,
456
+
t4_in: u64,
457
+
t5_in: u64,
458
+
t6_in: u64,
459
+
t7_in: u64,
460
+
t8_in: u64,
461
+
t9_in: u64,
462
+
t10: u64,
463
+
t11: u64,
464
+
t12: u64,
465
+
t13: u64,
466
+
t14: u64,
467
+
t15: u64,
468
+
t16: u64,
469
+
t17: u64,
470
+
t18: u64,
471
+
t19: u64,
472
+
) Fe26 {
473
+
var m2 = t0_in + t10 * 15632;
474
+
const r0 = m2 & field_base_mask;
475
+
m2 = (m2 >> field_base) + t1_in + t10 * 1024 + t11 * 15632;
476
+
const r1 = m2 & field_base_mask;
477
+
m2 = (m2 >> field_base) + t2_in + t11 * 1024 + t12 * 15632;
478
+
const r2 = m2 & field_base_mask;
479
+
m2 = (m2 >> field_base) + t3_in + t12 * 1024 + t13 * 15632;
480
+
const t3 = m2 & field_base_mask;
481
+
m2 = (m2 >> field_base) + t4_in + t13 * 1024 + t14 * 15632;
482
+
const t4 = m2 & field_base_mask;
483
+
m2 = (m2 >> field_base) + t5_in + t14 * 1024 + t15 * 15632;
484
+
const t5 = m2 & field_base_mask;
485
+
m2 = (m2 >> field_base) + t6_in + t15 * 1024 + t16 * 15632;
486
+
const t6 = m2 & field_base_mask;
487
+
m2 = (m2 >> field_base) + t7_in + t16 * 1024 + t17 * 15632;
488
+
const t7 = m2 & field_base_mask;
489
+
m2 = (m2 >> field_base) + t8_in + t17 * 1024 + t18 * 15632;
490
+
const t8 = m2 & field_base_mask;
491
+
m2 = (m2 >> field_base) + t9_in + t18 * 1024 + t19 * 68719492368;
492
+
const t9 = m2 & field_msb_mask;
493
+
m2 >>= field_msb_bits;
494
+
495
+
// final single-iteration fixup
496
+
const d0 = r0 + m2 * 977;
497
+
const d1 = (d0 >> field_base) + r1 + m2 * 64;
498
+
499
+
return .{ .n = .{
500
+
@truncate(d0 & field_base_mask),
501
+
@truncate(d1 & field_base_mask),
502
+
@truncate((d1 >> field_base) + r2),
503
+
@truncate(t3),
504
+
@truncate(t4),
505
+
@truncate(t5),
506
+
@truncate(t6),
507
+
@truncate(t7),
508
+
@truncate(t8),
509
+
@truncate(t9),
510
+
} };
511
+
}
512
+
513
+
/// Comptime conversion from integer literal to Fe26.
514
+
pub fn fromInt(comptime val: u256) Fe26 {
515
+
var b: [32]u8 = undefined;
516
+
std.mem.writeInt(u256, &b, val, .big);
517
+
return fromBytes(b, .big);
518
+
}
519
+
520
+
/// Modular inverse via Fermat's little theorem: a^(p-2) mod p.
521
+
/// Uses the dcrd addition chain: 255 squarings + 15 multiplications.
522
+
pub fn invert(self: Fe26) Fe26 {
523
+
const a = self;
524
+
const a2 = a.sq().mul(a); // a^(2^2 - 1)
525
+
const a3 = a2.sq().mul(a); // a^(2^3 - 1)
526
+
const a6 = sqn(a3, 3).mul(a3); // a^(2^6 - 1)
527
+
const a9 = sqn(a6, 3).mul(a3); // a^(2^9 - 1)
528
+
const a11 = sqn(a9, 2).mul(a2); // a^(2^11 - 1)
529
+
const a22 = sqn(a11, 11).mul(a11); // a^(2^22 - 1)
530
+
const a44 = sqn(a22, 22).mul(a22); // a^(2^44 - 1)
531
+
const a88 = sqn(a44, 44).mul(a44); // a^(2^88 - 1)
532
+
const a176 = sqn(a88, 88).mul(a88); // a^(2^176 - 1)
533
+
const a220 = sqn(a176, 44).mul(a44); // a^(2^220 - 1)
534
+
const a223 = sqn(a220, 3).mul(a3); // a^(2^223 - 1)
535
+
536
+
// final: a^(2^256 - 4294968275) = a^(p-2)
537
+
var r = sqn(a223, 23).mul(a22); // a^(2^246 - 4194305)
538
+
r = sqn(r, 5).mul(a); // a^(2^251 - 134217759)
539
+
r = sqn(r, 3).mul(a2); // a^(2^254 - 1073742069)
540
+
r = sqn(r, 2).mul(a); // a^(2^256 - 4294968275) = a^(p-2)
541
+
return r;
542
+
}
543
+
544
+
/// Helper: square n times.
545
+
fn sqn(f: Fe26, comptime n: usize) Fe26 {
546
+
var r = f;
547
+
for (0..n) |_| r = r.sq();
548
+
return r;
549
+
}
550
+
551
+
/// Convert from stdlib Fe to Fe26 via byte round-trip.
552
+
pub fn fromStdlib(fe: StdFe) Fe26 {
553
+
return fromBytes(fe.toBytes(.big), .big);
554
+
}
555
+
556
+
/// Convert to stdlib Fe via byte round-trip.
557
+
pub fn toStdlib(self: Fe26) StdFe {
558
+
return StdFe.fromBytes(self.normalize().toBytes(.big), .big) catch unreachable;
559
+
}
560
+
};
561
+
562
+
// ============================================================
563
+
// tests
564
+
// ============================================================
565
+
566
+
const testing = std.testing;
567
+
568
+
fn stdFeMul(a: StdFe, b: StdFe) StdFe {
569
+
return a.mul(b);
570
+
}
571
+
572
+
fn stdFeSq(a: StdFe) StdFe {
573
+
return a.sq();
574
+
}
575
+
576
+
test "fromBytes/toBytes round-trip" {
577
+
// generator point x-coordinate
578
+
const gx_bytes = Secp256k1.basePoint.affineCoordinates().x.toBytes(.big);
579
+
const fe = Fe26.fromBytes(gx_bytes, .big);
580
+
const out = fe.normalize().toBytes(.big);
581
+
try testing.expectEqualSlices(u8, &gx_bytes, &out);
582
+
}
583
+
584
+
test "fromBytes/toBytes little-endian round-trip" {
585
+
const gx_bytes = Secp256k1.basePoint.affineCoordinates().x.toBytes(.little);
586
+
const fe = Fe26.fromBytes(gx_bytes, .little);
587
+
const out = fe.normalize().toBytes(.little);
588
+
try testing.expectEqualSlices(u8, &gx_bytes, &out);
589
+
}
590
+
591
+
test "zero and one" {
592
+
try testing.expect(Fe26.zero.isZero());
593
+
try testing.expect(!Fe26.one.isZero());
594
+
try testing.expect(Fe26.one.equivalent(Fe26.one));
595
+
try testing.expect(!Fe26.zero.equivalent(Fe26.one));
596
+
}
597
+
598
+
test "add matches stdlib" {
599
+
const g = Secp256k1.basePoint.affineCoordinates();
600
+
const a = Fe26.fromStdlib(g.x);
601
+
const b = Fe26.fromStdlib(g.y);
602
+
const result = a.add(b).normalize().toBytes(.big);
603
+
const expected = g.x.add(g.y).toBytes(.big);
604
+
try testing.expectEqualSlices(u8, &expected, &result);
605
+
}
606
+
607
+
test "sub matches stdlib" {
608
+
const g = Secp256k1.basePoint.affineCoordinates();
609
+
const a = Fe26.fromStdlib(g.x);
610
+
const b = Fe26.fromStdlib(g.y);
611
+
const result = a.sub(b).normalize().toBytes(.big);
612
+
const expected = g.x.sub(g.y).toBytes(.big);
613
+
try testing.expectEqualSlices(u8, &expected, &result);
614
+
}
615
+
616
+
test "neg matches stdlib" {
617
+
const g = Secp256k1.basePoint.affineCoordinates();
618
+
const a = Fe26.fromStdlib(g.x);
619
+
const result = a.neg(1).normalize().toBytes(.big);
620
+
const expected = StdFe.zero.sub(g.x).toBytes(.big);
621
+
try testing.expectEqualSlices(u8, &expected, &result);
622
+
}
623
+
624
+
test "mul matches stdlib" {
625
+
const g = Secp256k1.basePoint.affineCoordinates();
626
+
const a = Fe26.fromStdlib(g.x);
627
+
const b = Fe26.fromStdlib(g.y);
628
+
const result = a.mul(b).normalize().toBytes(.big);
629
+
const expected = stdFeMul(g.x, g.y).toBytes(.big);
630
+
try testing.expectEqualSlices(u8, &expected, &result);
631
+
}
632
+
633
+
test "sq matches stdlib" {
634
+
const g = Secp256k1.basePoint.affineCoordinates();
635
+
const a = Fe26.fromStdlib(g.x);
636
+
const result = a.sq().normalize().toBytes(.big);
637
+
const expected = stdFeSq(g.x).toBytes(.big);
638
+
try testing.expectEqualSlices(u8, &expected, &result);
639
+
}
640
+
641
+
test "sq equals mul(a, a)" {
642
+
const g = Secp256k1.basePoint.affineCoordinates();
643
+
const a = Fe26.fromStdlib(g.x);
644
+
try testing.expect(a.sq().equivalent(a.mul(a)));
645
+
646
+
const b = Fe26.fromStdlib(g.y);
647
+
try testing.expect(b.sq().equivalent(b.mul(b)));
648
+
}
649
+
650
+
test "mul(a, invert(a)) == one" {
651
+
const g = Secp256k1.basePoint.affineCoordinates();
652
+
const a = Fe26.fromStdlib(g.x);
653
+
const inv = a.invert();
654
+
const product = a.mul(inv);
655
+
try testing.expect(product.equivalent(Fe26.one));
656
+
}
657
+
658
+
test "mul(a, invert(a)) == one for y" {
659
+
const g = Secp256k1.basePoint.affineCoordinates();
660
+
const a = Fe26.fromStdlib(g.y);
661
+
const inv = a.invert();
662
+
const product = a.mul(inv);
663
+
try testing.expect(product.equivalent(Fe26.one));
664
+
}
665
+
666
+
test "edge: p-1" {
667
+
// p - 1 = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2E
668
+
const p_minus_1 = Fe26.fromInt(
669
+
0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2E,
670
+
);
671
+
try testing.expect(!p_minus_1.isZero());
672
+
// p-1 + 1 should normalize to 0
673
+
try testing.expect(p_minus_1.add(Fe26.one).isZero());
674
+
}
675
+
676
+
test "fromInt constants" {
677
+
// beta (endomorphism constant) as Fe26
678
+
const beta = Fe26.fromInt(
679
+
55594575648329892869085402983802832744385952214688224221778511981742606582254,
680
+
);
681
+
try testing.expect(!beta.isZero());
682
+
683
+
// verify beta^3 == 1 (mod p) — cube root of unity
684
+
const beta_cubed = beta.sq().mul(beta);
685
+
try testing.expect(beta_cubed.equivalent(Fe26.one));
686
+
}
687
+
688
+
test "dbl matches add(self)" {
689
+
const g = Secp256k1.basePoint.affineCoordinates();
690
+
const a = Fe26.fromStdlib(g.x);
691
+
try testing.expect(a.dbl().equivalent(a.add(a)));
692
+
}
693
+
694
+
test "mulInt matches repeated add" {
695
+
const g = Secp256k1.basePoint.affineCoordinates();
696
+
const a = Fe26.fromStdlib(g.x);
697
+
try testing.expect(a.mulInt(3).equivalent(a.add(a).add(a)));
698
+
}
699
+
700
+
test "multiple mul/sq operations stay consistent" {
701
+
// compute (gx^2 * gy) via Fe26 and stdlib, compare
702
+
const g = Secp256k1.basePoint.affineCoordinates();
703
+
const gx26 = Fe26.fromStdlib(g.x);
704
+
const gy26 = Fe26.fromStdlib(g.y);
705
+
706
+
const result = gx26.sq().mul(gy26);
707
+
const expected = stdFeSq(g.x).mul(g.y);
708
+
try testing.expectEqualSlices(u8, &expected.toBytes(.big), &result.normalize().toBytes(.big));
709
+
}
710
+
711
+
test "normalize idempotent" {
712
+
const g = Secp256k1.basePoint.affineCoordinates();
713
+
const a = Fe26.fromStdlib(g.x);
714
+
const n1 = a.normalize();
715
+
const n2 = n1.normalize();
716
+
try testing.expectEqualSlices(u8, &std.mem.toBytes(n1.n), &std.mem.toBytes(n2.n));
717
+
}
+211
-59
src/jacobian.zig
+211
-59
src/jacobian.zig
···
1
1
const std = @import("std");
2
2
const Secp256k1 = std.crypto.ecc.Secp256k1;
3
-
const Fe = Secp256k1.Fe;
4
-
const AffineCoordinates = @TypeOf(Secp256k1.basePoint.affineCoordinates());
3
+
const Fe26 = @import("field.zig").Fe26;
4
+
const StdFe = Secp256k1.Fe;
5
+
const StdAffineCoordinates = @TypeOf(Secp256k1.basePoint.affineCoordinates());
6
+
7
+
/// Affine point in Fe26 representation.
8
+
pub const AffinePoint26 = struct {
9
+
x: Fe26,
10
+
y: Fe26,
11
+
12
+
pub const identity: AffinePoint26 = .{ .x = Fe26.zero, .y = Fe26.zero };
5
13
6
-
/// Jacobian point on secp256k1: affine (X/Z², Y/Z³).
14
+
pub fn neg(self: AffinePoint26) AffinePoint26 {
15
+
return .{ .x = self.x, .y = self.y.neg(1) };
16
+
}
17
+
18
+
pub fn fromStdlib(ac: StdAffineCoordinates) AffinePoint26 {
19
+
return .{
20
+
.x = Fe26.fromBytes(ac.x.toBytes(.big), .big),
21
+
.y = Fe26.fromBytes(ac.y.toBytes(.big), .big),
22
+
};
23
+
}
24
+
25
+
pub fn fromStdlibIdentity(ac: StdAffineCoordinates) AffinePoint26 {
26
+
// check if this is the identity element (x=0, y=0 in stdlib)
27
+
if (ac.x.isZero() and ac.y.isZero()) return identity;
28
+
return fromStdlib(ac);
29
+
}
30
+
};
31
+
32
+
/// Jacobian point on secp256k1 using Fe26: affine (X/Z², Y/Z³).
7
33
/// Uses a=0 specialized formulas for fewer field operations than
8
34
/// the stdlib's complete projective formulas.
9
35
pub const JacobianPoint = struct {
10
-
x: Fe,
11
-
y: Fe,
12
-
z: Fe,
36
+
x: Fe26,
37
+
y: Fe26,
38
+
z: Fe26,
13
39
14
-
pub const identity: JacobianPoint = .{ .x = Fe.zero, .y = Fe.one, .z = Fe.zero };
40
+
pub const identity: JacobianPoint = .{ .x = Fe26.zero, .y = Fe26.one, .z = Fe26.zero };
15
41
16
-
/// Convert affine (x, y) to Jacobian (x, y, 1).
17
-
pub fn fromAffine(p: AffineCoordinates) JacobianPoint {
18
-
return .{ .x = p.x, .y = p.y, .z = Fe.one };
42
+
/// Convert Fe26 affine (x, y) to Jacobian (x, y, 1).
43
+
pub fn fromAffine(p: AffinePoint26) JacobianPoint {
44
+
return .{ .x = p.x, .y = p.y, .z = Fe26.one };
19
45
}
20
46
21
-
/// Convert Jacobian to projective (stdlib Secp256k1 point).
22
-
/// Jacobian: x = X/Z², y = Y/Z³. Projective: x = X'/Z', y = Y'/Z'.
23
-
/// Set X' = X*Z, Y' = Y, Z' = Z³.
24
-
pub fn toProjective(self: JacobianPoint) Secp256k1 {
25
-
if (self.z.isZero()) return Secp256k1.identityElement;
26
-
return .{
27
-
.x = self.x.mul(self.z),
28
-
.y = self.y,
29
-
.z = self.z.sq().mul(self.z),
30
-
};
47
+
/// Convert stdlib affine to Jacobian via Fe26.
48
+
pub fn fromStdlibAffine(ac: StdAffineCoordinates) JacobianPoint {
49
+
return fromAffine(AffinePoint26.fromStdlib(ac));
31
50
}
32
51
33
-
/// Doubling for a=0 curves (secp256k1). 2M + 5S.
52
+
/// Doubling for a=0 curves (secp256k1). 2M + 5S + 3 normalize.
34
53
/// EFD dbl-2009-l with a=0 specialization.
35
54
pub fn dbl(self: JacobianPoint) JacobianPoint {
36
55
if (self.z.isZero()) return self;
37
56
38
-
const a = self.x.sq(); // X1²
39
-
const b = self.y.sq(); // Y1²
40
-
const c = b.sq(); // Y1⁴
41
-
const d = self.x.add(b).sq().sub(a).sub(c).dbl(); // 4*X1*Y1²
42
-
const e = a.dbl().add(a); // 3*X1²
43
-
const f = e.sq(); // (3*X1²)²
44
-
const x3 = f.sub(d.dbl()); // F - 2*D
45
-
const c8 = c.dbl().dbl().dbl(); // 8*Y1⁴
46
-
const y3 = e.mul(d.sub(x3)).sub(c8); // E*(D-X3) - 8*C
47
-
const z3 = self.y.mul(self.z).dbl(); // 2*Y1*Z1
57
+
const a = self.x.sq(); // mag 1
58
+
const b = self.y.sq(); // mag 1
59
+
const c = b.sq(); // mag 1
60
+
const xb = self.x.add(b); // mag 2
61
+
const d = xb.sq().subMag(a, 1).subMag(c, 1).dbl().normalize();
62
+
const e = a.dbl().add(a); // 3a = mag 3
63
+
const f = e.sq();
64
+
const x3 = f.subMag(d.dbl(), 2).normalize();
65
+
const c8 = c.dbl().dbl().dbl(); // mag 8
66
+
const y3 = e.mul(d.subMag(x3, 1)).subMag(c8, 8).normalize();
67
+
const z3 = self.y.mul(self.z).dbl(); // mag 2
48
68
49
69
return .{ .x = x3, .y = y3, .z = z3 };
50
70
}
51
71
52
-
/// Mixed addition: Jacobian + Affine → Jacobian. 7M + 4S.
53
-
/// EFD madd-2007-bl.
54
-
pub fn addMixed(self: JacobianPoint, q: AffineCoordinates) JacobianPoint {
72
+
/// Mixed addition: Jacobian + Affine → Jacobian. 7M + 4S + 4 normalize.
73
+
/// EFD madd-2007-bl. Falls back to dbl() when both points are equal.
74
+
pub fn addMixed(self: JacobianPoint, q: AffinePoint26) JacobianPoint {
55
75
if (self.z.isZero()) return fromAffine(q);
56
76
57
-
const z1z1 = self.z.sq(); // Z1²
58
-
const qx_z2 = q.x.mul(z1z1); // x2*Z1²
59
-
const s2 = q.y.mul(self.z.mul(z1z1)); // y2*Z1³
60
-
const h = qx_z2.sub(self.x); // U2 - X1
61
-
const hh = h.sq(); // H²
62
-
const i = hh.dbl().dbl(); // 4*H²
63
-
const j = h.mul(i); // H*4*H²
64
-
const r = s2.sub(self.y).dbl(); // 2*(S2 - Y1)
65
-
const v = self.x.mul(i); // X1*4*H²
66
-
const x3 = r.sq().sub(j).sub(v.dbl()); // r²-J-2*V
67
-
const y3 = r.mul(v.sub(x3)).sub(self.y.mul(j).dbl()); // r*(V-X3)-2*Y1*J
68
-
const z3 = self.z.add(h).sq().sub(z1z1).sub(hh); // (Z1+H)²-Z1²-H²
77
+
const z1z1 = self.z.sq(); // mag 1
78
+
const qx_z2 = q.x.mul(z1z1); // mag 1
79
+
const s2 = q.y.mul(self.z.mul(z1z1)); // mag 1
80
+
const h = qx_z2.subMag(self.x, 1).normalize();
81
+
82
+
// Degenerate case: same x-coordinate
83
+
if (h.isZero()) {
84
+
if (s2.subMag(self.y, 1).normalize().isZero()) return self.dbl();
85
+
return identity;
86
+
}
87
+
88
+
const hh = h.sq();
89
+
const i = hh.dbl().dbl(); // mag 4
90
+
const j = h.mul(i);
91
+
const r = s2.subMag(self.y, 1).dbl().normalize();
92
+
const v = self.x.mul(i);
93
+
const x3 = r.sq().subMag(j, 1).subMag(v.dbl(), 2).normalize();
94
+
const y3 = r.mul(v.subMag(x3, 1)).subMag(self.y.mul(j).dbl(), 2).normalize();
95
+
const z3 = self.z.add(h).sq().subMag(z1z1, 1).subMag(hh, 1);
69
96
70
97
return .{ .x = x3, .y = y3, .z = z3 };
71
98
}
72
99
73
100
/// Mixed subtraction: Jacobian - Affine → Jacobian.
74
-
pub fn subMixed(self: JacobianPoint, q: AffineCoordinates) JacobianPoint {
101
+
pub fn subMixed(self: JacobianPoint, q: AffinePoint26) JacobianPoint {
75
102
return self.addMixed(q.neg());
76
103
}
104
+
105
+
/// Negate the Y coordinate: returns (X, -Y, Z).
106
+
pub fn negY(self: JacobianPoint) JacobianPoint {
107
+
return .{ .x = self.x, .y = self.y.normalize().neg(1), .z = self.z };
108
+
}
109
+
110
+
/// Full Jacobian addition: J + J → J. ~12M + 4S.
111
+
/// Needed for combining r1 + r2 at end of verify.
112
+
pub fn add(self: JacobianPoint, other: JacobianPoint) JacobianPoint {
113
+
if (self.z.isZero()) return other;
114
+
if (other.z.isZero()) return self;
115
+
116
+
const z1z1 = self.z.sq();
117
+
const z2z2 = other.z.sq();
118
+
const p1 = self.x.mul(z2z2); // X1*Z2²
119
+
const p2 = other.x.mul(z1z1); // X2*Z1²
120
+
const s1 = self.y.mul(other.z.mul(z2z2)); // Y1*Z2³
121
+
const s2 = other.y.mul(self.z.mul(z1z1)); // Y2*Z1³
122
+
123
+
// h = p2 - p1 ; both mag 1
124
+
const h = p2.subMag(p1, 1).normalize();
125
+
const r = s2.subMag(s1, 1).normalize();
126
+
127
+
// check for doubling case (h == 0 && r == 0)
128
+
if (h.isZero()) {
129
+
if (r.isZero()) return self.dbl();
130
+
return identity;
131
+
}
132
+
133
+
const hh = h.sq(); // mag 1
134
+
const hhh = h.mul(hh); // mag 1
135
+
const v = p1.mul(hh); // mag 1
136
+
const r2 = r.sq(); // mag 1
137
+
// x3 = r² - hhh - 2v ; sub(1)=3, sub(2v=2, 2)=6
138
+
const x3 = r2.subMag(hhh, 1).subMag(v.dbl(), 2).normalize();
139
+
// y3 = r*(v-x3) - s1*hhh ; v.sub(x3,1)→3, mul=1; s1.mul(hhh)=1; sub(1)=3
140
+
const y3 = r.mul(v.subMag(x3, 1)).subMag(s1.mul(hhh), 1).normalize();
141
+
const z3 = self.z.mul(other.z).mul(h); // mag 1
142
+
143
+
return .{ .x = x3, .y = y3, .z = z3 };
144
+
}
145
+
146
+
/// Compare signature r against R's x-coordinate in Jacobian space.
147
+
/// Checks r * Z² == X (mod p), avoiding field inversion.
148
+
/// Returns true if the x-coordinates match.
149
+
pub fn jacobianCompare(self: JacobianPoint, r_fe: Fe26) bool {
150
+
// Jacobian: x_affine = X / Z²
151
+
// So check: r * Z² == X
152
+
const z2 = self.z.sq();
153
+
const rz2 = r_fe.mul(z2);
154
+
if (rz2.equivalent(self.x)) return true;
155
+
156
+
// rare overflow case: x_affine could be in [n, p)
157
+
// check (r + n) * Z² == X
158
+
const n_fe26 = comptime Fe26.fromInt(Secp256k1.scalar.field_order);
159
+
const p_minus_n: u256 = StdFe.field_order - Secp256k1.scalar.field_order;
160
+
const r_norm = r_fe.normalize();
161
+
const r_bytes = r_norm.toBytes(.big);
162
+
const r_int = std.mem.readInt(u256, &r_bytes, .big);
163
+
if (r_int < p_minus_n) {
164
+
const rn_z2 = r_fe.add(n_fe26).mul(z2);
165
+
if (rn_z2.equivalent(self.x)) return true;
166
+
}
167
+
168
+
return false;
169
+
}
170
+
171
+
/// Convert to stdlib projective for backward compatibility.
172
+
/// Jacobian: x = X/Z², y = Y/Z³. Projective: x = X'/Z', y = Y'/Z'.
173
+
/// Set X' = X*Z, Y' = Y, Z' = Z³.
174
+
pub fn toProjective(self: JacobianPoint) Secp256k1 {
175
+
if (self.z.isZero()) return Secp256k1.identityElement;
176
+
return .{
177
+
.x = self.x.mul(self.z).toStdlib(),
178
+
.y = self.y.toStdlib(),
179
+
.z = self.z.sq().mul(self.z).toStdlib(),
180
+
};
181
+
}
77
182
};
78
183
184
+
// ============================================================
185
+
// tests — verify Fe26 Jacobian matches stdlib results
186
+
// ============================================================
187
+
188
+
fn toStdlibAffine(j: JacobianPoint) StdAffineCoordinates {
189
+
return j.toProjective().affineCoordinates();
190
+
}
191
+
79
192
test "dbl matches stdlib" {
80
193
const G = Secp256k1.basePoint;
81
194
const g_affine = G.affineCoordinates();
82
195
83
-
const j_2g = JacobianPoint.fromAffine(g_affine).dbl();
84
-
const j_2g_affine = j_2g.toProjective().affineCoordinates();
196
+
const j_2g = JacobianPoint.fromStdlibAffine(g_affine).dbl();
197
+
const j_2g_affine = toStdlibAffine(j_2g);
85
198
86
199
const stdlib_2g_affine = G.dbl().affineCoordinates();
87
200
···
92
205
test "addMixed matches stdlib add" {
93
206
const G = Secp256k1.basePoint;
94
207
const g_affine = G.affineCoordinates();
208
+
const g26 = AffinePoint26.fromStdlib(g_affine);
95
209
96
210
// 2G + G = 3G
97
-
const j_3g = JacobianPoint.fromAffine(g_affine).dbl().addMixed(g_affine);
98
-
const j_3g_affine = j_3g.toProjective().affineCoordinates();
211
+
const j_3g = JacobianPoint.fromAffine(g26).dbl().addMixed(g26);
212
+
const j_3g_affine = toStdlibAffine(j_3g);
99
213
100
214
const stdlib_3g_affine = G.dbl().add(G).affineCoordinates();
101
215
···
106
220
test "identity + affine = affine" {
107
221
const G = Secp256k1.basePoint;
108
222
const g_affine = G.affineCoordinates();
223
+
const g26 = AffinePoint26.fromStdlib(g_affine);
109
224
110
-
const result = JacobianPoint.identity.addMixed(g_affine);
111
-
const result_affine = result.toProjective().affineCoordinates();
225
+
const result = JacobianPoint.identity.addMixed(g26);
226
+
const result_affine = toStdlibAffine(result);
112
227
113
228
try std.testing.expect(result_affine.x.equivalent(g_affine.x));
114
229
try std.testing.expect(result_affine.y.equivalent(g_affine.y));
···
117
232
test "subMixed is inverse of addMixed" {
118
233
const G = Secp256k1.basePoint;
119
234
const g_affine = G.affineCoordinates();
235
+
const g26 = AffinePoint26.fromStdlib(g_affine);
120
236
121
237
// 3G - G = 2G
122
-
const j_3g = JacobianPoint.fromAffine(g_affine).dbl().addMixed(g_affine);
123
-
const j_2g = j_3g.subMixed(g_affine);
124
-
const j_2g_affine = j_2g.toProjective().affineCoordinates();
238
+
const j_3g = JacobianPoint.fromAffine(g26).dbl().addMixed(g26);
239
+
const j_2g = j_3g.subMixed(g26);
240
+
const j_2g_affine = toStdlibAffine(j_2g);
125
241
126
242
const stdlib_2g_affine = G.dbl().affineCoordinates();
127
243
···
134
250
const g_affine = G.affineCoordinates();
135
251
136
252
// 16G via 4 doublings
137
-
var j = JacobianPoint.fromAffine(g_affine);
138
-
j = j.dbl().dbl().dbl().dbl();
139
-
const j_affine = j.toProjective().affineCoordinates();
253
+
const j = JacobianPoint.fromStdlibAffine(g_affine).dbl().dbl().dbl().dbl();
254
+
const j_affine = toStdlibAffine(j);
140
255
141
256
const stdlib_affine = G.dbl().dbl().dbl().dbl().affineCoordinates();
142
257
143
258
try std.testing.expect(j_affine.x.equivalent(stdlib_affine.x));
144
259
try std.testing.expect(j_affine.y.equivalent(stdlib_affine.y));
145
260
}
261
+
262
+
test "full Jacobian add matches stdlib" {
263
+
const G = Secp256k1.basePoint;
264
+
const g_affine = G.affineCoordinates();
265
+
const g26 = AffinePoint26.fromStdlib(g_affine);
266
+
267
+
// 2G + 3G = 5G
268
+
const j_2g = JacobianPoint.fromAffine(g26).dbl();
269
+
const j_3g = j_2g.addMixed(g26);
270
+
const j_5g = j_2g.add(j_3g);
271
+
const j_5g_affine = toStdlibAffine(j_5g);
272
+
273
+
const stdlib_5g = G.dbl().add(G.dbl().add(G));
274
+
const stdlib_5g_affine = stdlib_5g.affineCoordinates();
275
+
276
+
try std.testing.expect(j_5g_affine.x.equivalent(stdlib_5g_affine.x));
277
+
try std.testing.expect(j_5g_affine.y.equivalent(stdlib_5g_affine.y));
278
+
}
279
+
280
+
test "add with identity" {
281
+
const G = Secp256k1.basePoint;
282
+
const g_affine = G.affineCoordinates();
283
+
const g26 = AffinePoint26.fromStdlib(g_affine);
284
+
285
+
const j_g = JacobianPoint.fromAffine(g26);
286
+
const result = j_g.add(JacobianPoint.identity);
287
+
const result_affine = toStdlibAffine(result);
288
+
289
+
try std.testing.expect(result_affine.x.equivalent(g_affine.x));
290
+
try std.testing.expect(result_affine.y.equivalent(g_affine.y));
291
+
292
+
const result2 = JacobianPoint.identity.add(j_g);
293
+
const result2_affine = toStdlibAffine(result2);
294
+
295
+
try std.testing.expect(result2_affine.x.equivalent(g_affine.x));
296
+
try std.testing.expect(result2_affine.y.equivalent(g_affine.y));
297
+
}
+35
-23
src/point.zig
+35
-23
src/point.zig
···
1
1
const std = @import("std");
2
-
const Secp256k1 = std.crypto.ecc.Secp256k1;
2
+
const jacobian_mod = @import("jacobian.zig");
3
+
const JacobianPoint = jacobian_mod.JacobianPoint;
4
+
const AffinePoint26 = jacobian_mod.AffinePoint26;
3
5
const endo = @import("endo.zig");
6
+
const affine_mod = @import("affine.zig");
7
+
const Secp256k1 = std.crypto.ecc.Secp256k1;
4
8
5
9
/// Build a precomputation table: [identity, p, 2p, 3p, ..., count*p].
6
-
/// Matches stdlib's algorithm (double-or-add construction).
7
-
pub fn precompute(p: Secp256k1, comptime count: usize) [1 + count]Secp256k1 {
8
-
var pc: [1 + count]Secp256k1 = undefined;
9
-
pc[0] = Secp256k1.identityElement;
10
-
pc[1] = p;
10
+
/// Uses Jacobian arithmetic for efficient computation.
11
+
pub fn precompute(p: AffinePoint26, comptime count: usize) [1 + count]JacobianPoint {
12
+
var pc: [1 + count]JacobianPoint = undefined;
13
+
pc[0] = JacobianPoint.identity;
14
+
pc[1] = JacobianPoint.fromAffine(p);
11
15
var i: usize = 2;
12
16
while (i <= count) : (i += 1) {
13
-
pc[i] = if (i % 2 == 0) pc[i / 2].dbl() else pc[i - 1].add(p);
17
+
pc[i] = if (i % 2 == 0) pc[i / 2].dbl() else pc[i - 1].addMixed(p);
14
18
}
15
19
return pc;
16
20
}
17
21
18
-
/// Apply the endomorphism to each entry in a precompute table.
22
+
/// Apply the endomorphism to each entry in an affine precompute table.
19
23
/// phi(n*P) = n*phi(P), so transforming the table gives a valid table for phi(P).
20
-
pub fn phiTable(comptime n: usize, table: [n]Secp256k1) [n]Secp256k1 {
21
-
var result: [n]Secp256k1 = undefined;
24
+
pub fn phiTableAffine(comptime n: usize, table: [n]AffinePoint26) [n]AffinePoint26 {
25
+
var result: [n]AffinePoint26 = undefined;
22
26
for (0..n) |i| {
23
-
result[i] = endo.phi(table[i]);
27
+
result[i] = endo.phiAffine(table[i]);
24
28
}
25
29
return result;
26
30
}
27
31
28
-
const AffineCoordinates = @TypeOf(Secp256k1.basePoint.affineCoordinates());
29
-
30
-
/// Apply the endomorphism to each entry in an affine precompute table.
31
-
pub fn phiTableAffine(comptime n: usize, table: [n]AffineCoordinates) [n]AffineCoordinates {
32
-
var result: [n]AffineCoordinates = undefined;
32
+
/// Apply the endomorphism to each entry in a Jacobian precompute table.
33
+
pub fn phiTableJacobian(comptime n: usize, table: [n]JacobianPoint) [n]JacobianPoint {
34
+
var result: [n]JacobianPoint = undefined;
33
35
for (0..n) |i| {
34
-
result[i] = endo.phiAffine(table[i]);
36
+
result[i] = endo.phiJacobian(table[i]);
35
37
}
36
38
return result;
37
39
}
···
57
59
58
60
test "precompute table correctness" {
59
61
const G = Secp256k1.basePoint;
60
-
const table = precompute(G, 8);
62
+
const g_affine = G.affineCoordinates();
63
+
const g26 = AffinePoint26.fromStdlib(g_affine);
64
+
65
+
const table = precompute(g26, 8);
61
66
62
67
// table[0] should be identity
63
-
try std.testing.expect(table[0].x.equivalent(Secp256k1.identityElement.x));
64
68
try std.testing.expect(table[0].z.isZero());
65
69
66
-
// table[1] should be G
67
-
try std.testing.expect(Secp256k1.equivalent(table[1], G));
70
+
// table[1] should be G — compare via toProjective
71
+
const t1_affine = table[1].toProjective().affineCoordinates();
72
+
try std.testing.expect(t1_affine.x.equivalent(g_affine.x));
73
+
try std.testing.expect(t1_affine.y.equivalent(g_affine.y));
68
74
69
75
// table[2] should be 2G
70
-
try std.testing.expect(Secp256k1.equivalent(table[2], G.dbl()));
76
+
const t2_affine = table[2].toProjective().affineCoordinates();
77
+
const stdlib_2g = G.dbl().affineCoordinates();
78
+
try std.testing.expect(t2_affine.x.equivalent(stdlib_2g.x));
79
+
try std.testing.expect(t2_affine.y.equivalent(stdlib_2g.y));
71
80
72
81
// table[3] should be 3G = 2G + G
73
-
try std.testing.expect(Secp256k1.equivalent(table[3], G.dbl().add(G)));
82
+
const t3_affine = table[3].toProjective().affineCoordinates();
83
+
const stdlib_3g = G.dbl().add(G).affineCoordinates();
84
+
try std.testing.expect(t3_affine.x.equivalent(stdlib_3g.x));
85
+
try std.testing.expect(t3_affine.y.equivalent(stdlib_3g.y));
74
86
}
75
87
76
88
test "slide produces valid digits" {
+2
src/root.zig
+2
src/root.zig
···
7
7
const verify_mod = @import("verify.zig");
8
8
9
9
// re-export internals for testing
10
+
pub const field = @import("field.zig");
10
11
pub const endo = @import("endo.zig");
11
12
pub const point = @import("point.zig");
12
13
pub const verify = @import("verify.zig");
···
123
124
};
124
125
125
126
test {
127
+
_ = @import("field.zig");
126
128
_ = @import("endo.zig");
127
129
_ = @import("point.zig");
128
130
_ = @import("verify.zig");
+59
-46
src/verify.zig
+59
-46
src/verify.zig
···
1
1
const std = @import("std");
2
-
const mem = std.mem;
3
2
const Secp256k1 = std.crypto.ecc.Secp256k1;
4
-
const Fe = Secp256k1.Fe;
5
-
const AffineCoordinates = @TypeOf(Secp256k1.basePoint.affineCoordinates());
3
+
const Fe26 = @import("field.zig").Fe26;
6
4
const scalar = Secp256k1.scalar;
7
5
8
6
const endo = @import("endo.zig");
···
10
8
const affine_mod = @import("affine.zig");
11
9
const jacobian_mod = @import("jacobian.zig");
12
10
const JacobianPoint = jacobian_mod.JacobianPoint;
11
+
const AffinePoint26 = jacobian_mod.AffinePoint26;
13
12
14
13
pub const VerifyError = std.crypto.errors.IdentityElementError ||
15
14
std.crypto.errors.NonCanonicalError ||
16
15
error{SignatureVerificationFailed};
17
16
18
-
/// Precomputed base point table: G_TABLE[i][j] = j * 256^i * G in affine.
19
-
/// 16 subtables × 256 entries = 4096 affine points (~256KB).
17
+
/// Precomputed base point table: G_TABLE[i][j] = j * 256^i * G in Fe26 affine.
18
+
/// 16 subtables × 256 entries = 4096 affine points.
20
19
/// Enables u1*G computation via byte-at-a-time lookup with zero doublings.
21
-
const G_TABLE: [16][256]AffineCoordinates = blk: {
22
-
@setEvalBranchQuota(50_000_000);
20
+
const G_TABLE: [16][256]AffinePoint26 = blk: {
21
+
@setEvalBranchQuota(100_000_000);
23
22
break :blk affine_mod.buildByteTable(Secp256k1.basePoint);
24
23
};
25
24
26
-
/// Curve order as a field element, for the projective overflow check.
27
-
const n_fe = Fe.fromInt(scalar.field_order) catch unreachable;
28
-
29
-
/// p - n: if r (as integer) < this, then r + n < p and overflow is possible.
30
-
const p_minus_n: u256 = Fe.field_order - scalar.field_order;
31
-
32
25
/// Reduce a 32-byte big-endian value to a scalar (mod n).
33
26
fn reduceToScalar(h: [32]u8) scalar.Scalar {
34
27
var xs = [_]u8{0} ** 48;
···
38
31
39
32
/// Verify an ECDSA signature using precomputed tables and Jacobian arithmetic.
40
33
///
41
-
/// Optimizations over v0.1:
34
+
/// Optimizations:
42
35
/// 1. u1*G via precomputed byte table: ~32 mixed adds, zero doublings
43
36
/// 2. u2*Q via 2-way Jacobian Shamir: cheaper dbl (2M+5S) and mixed add (7M+4S)
44
-
/// 3. Endomorphism + projective comparison (carried from v0.1)
37
+
/// 3. Endomorphism + Jacobian comparison (no field inversion)
38
+
/// 4. All arithmetic in Fe26 (10×26-bit) — no Montgomery form overhead
45
39
pub fn verify(sig_r: [32]u8, sig_s: [32]u8, msg_hash: [32]u8, public_key: Secp256k1) VerifyError!void {
46
40
// parse and validate r, s
47
41
const r_sc = scalar.Scalar.fromBytes(sig_r, .big) catch return error.SignatureVerificationFailed;
···
87
81
const r1 = basePointMul(split_u1.r1, neg_g, split_u1.r2, neg_g_phi);
88
82
89
83
// 2. u2*Q via 2-way Jacobian Shamir
90
-
const pk_pc = point.precompute(public_key, 8);
91
-
const pk_affine = affine_mod.batchToAffine(9, pk_pc);
84
+
const pk_affine26 = AffinePoint26.fromStdlib(public_key.affineCoordinates());
85
+
const pk_jac = point.precompute(pk_affine26, 8);
86
+
const pk_affine = affine_mod.batchToAffine(9, pk_jac);
92
87
const pk_phi_affine = point.phiTableAffine(9, pk_affine);
93
88
const r2 = publicKeyMul(split_u2, neg_p, neg_p_phi, pk_affine, pk_phi_affine);
94
89
95
-
// 3. combine results: convert Jacobian → projective, add via stdlib
96
-
const p1 = r1.toProjective();
97
-
const p2 = r2.toProjective();
98
-
const q = p1.add(p2);
90
+
// 3. combine results in Jacobian, compare without inversion
91
+
const q = r1.add(r2);
99
92
100
-
// reject identity (point at infinity has no valid x-coordinate)
101
-
q.rejectIdentity() catch return error.SignatureVerificationFailed;
93
+
// reject identity
94
+
if (q.z.isZero()) return error.SignatureVerificationFailed;
102
95
103
-
// projective comparison: check x(R) mod n == r
104
-
if (!projectiveCompare(sig_r, q)) {
96
+
// Jacobian comparison: r * Z² == X (mod p)
97
+
const r_fe = Fe26.fromBytes(sig_r, .big);
98
+
if (!q.jacobianCompare(r_fe)) {
105
99
return error.SignatureVerificationFailed;
106
100
}
107
101
}
···
138
132
split: endo.SplitScalar,
139
133
neg_p: bool,
140
134
neg_p_phi: bool,
141
-
pk_affine: [9]AffineCoordinates,
142
-
pk_phi_affine: [9]AffineCoordinates,
135
+
pk_affine: [9]AffinePoint26,
136
+
pk_phi_affine: [9]AffinePoint26,
143
137
) JacobianPoint {
144
138
const e1 = point.slide(split.r1);
145
139
const e2 = point.slide(split.r2);
···
157
151
158
152
/// Add/subtract a table entry based on signed digit and negation flag.
159
153
/// Variable-time: branches on digit value (safe for public verification).
160
-
inline fn addSlot(q: JacobianPoint, table: *const [9]AffineCoordinates, slot: i8, negate: bool) JacobianPoint {
154
+
inline fn addSlot(q: JacobianPoint, table: *const [9]AffinePoint26, slot: i8, negate: bool) JacobianPoint {
161
155
var s = slot;
162
156
if (negate) s = -s;
163
157
if (s > 0) {
···
168
162
return q;
169
163
}
170
164
171
-
/// Compare signature r against R's x-coordinate in projective space.
172
-
/// The stdlib's Secp256k1 uses projective coordinates: x_affine = X / Z.
173
-
/// This avoids the expensive field inversion (~240 muls) by checking r * Z == X (mod p).
174
-
fn projectiveCompare(r_bytes: [32]u8, result: Secp256k1) bool {
175
-
// r < n < p, so this always succeeds
176
-
const r_fe = Fe.fromBytes(r_bytes, .big) catch return false;
177
-
const rz = r_fe.mul(result.z);
165
+
test "G_TABLE[0][1] is the base point" {
166
+
const G = Secp256k1.basePoint;
167
+
const g_affine = G.affineCoordinates();
178
168
179
-
// common case: x_affine < n, so r * Z == X (mod p)
180
-
if (result.x.equivalent(rz)) return true;
169
+
// G_TABLE[0][1] should be 1*G
170
+
const table_g = G_TABLE[0][1];
171
+
const table_g_x = table_g.x.toStdlib();
172
+
const table_g_y = table_g.y.toStdlib();
181
173
182
-
// rare case (probability ~2^-128): x_affine in [n, p)
183
-
// only possible if r + n < p
184
-
const r_int = mem.readInt(u256, &r_bytes, .big);
185
-
if (r_int < p_minus_n) {
186
-
const rn_z = r_fe.add(n_fe).mul(result.z);
187
-
if (result.x.equivalent(rn_z)) return true;
188
-
}
174
+
try std.testing.expect(table_g_x.equivalent(g_affine.x));
175
+
try std.testing.expect(table_g_y.equivalent(g_affine.y));
176
+
}
189
177
190
-
return false;
178
+
test "G_TABLE[0][2] is 2G" {
179
+
const G = Secp256k1.basePoint;
180
+
const g2 = G.dbl().affineCoordinates();
181
+
182
+
const table_2g = G_TABLE[0][2];
183
+
const table_2g_x = table_2g.x.toStdlib();
184
+
const table_2g_y = table_2g.y.toStdlib();
185
+
186
+
try std.testing.expect(table_2g_x.equivalent(g2.x));
187
+
try std.testing.expect(table_2g_y.equivalent(g2.y));
188
+
}
189
+
190
+
test "comptime addMixed gives correct 2G" {
191
+
@setEvalBranchQuota(100_000_000);
192
+
const g_affine = comptime AffinePoint26.fromStdlib(Secp256k1.basePoint.affineCoordinates());
193
+
const j_g = comptime JacobianPoint.fromAffine(g_affine);
194
+
const j_2g = comptime j_g.addMixed(g_affine);
195
+
196
+
// convert to affine via toProjective
197
+
const proj = j_2g.toProjective();
198
+
const affine = proj.affineCoordinates();
199
+
200
+
const expected = Secp256k1.basePoint.dbl().affineCoordinates();
201
+
202
+
try std.testing.expect(affine.x.equivalent(expected.x));
203
+
try std.testing.expect(affine.y.equivalent(expected.y));
191
204
}
192
205
193
206
test "projective x-coordinate: X == x_affine * Z" {
194
207
// verify the coordinate convention: x_affine = X / Z (projective, not Jacobian)
195
208
const s = comptime blk: {
196
209
var buf: [32]u8 = undefined;
197
-
mem.writeInt(u256, &buf, 12345, .little);
210
+
std.mem.writeInt(u256, &buf, 12345, .little);
198
211
break :blk buf;
199
212
};
200
213
const P = try Secp256k1.basePoint.mul(s, .little);