feat: 5×52-bit field, Fermat scalar inversion, type cleanup

+14 -14

src/affine.zig

··· 1 const std = @import("std"); 2 const Secp256k1 = std.crypto.ecc.Secp256k1; 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; 7 8 /// Batch-convert Jacobian points to affine using Montgomery's trick. 9 /// Uses 1 field inversion + 3*(n-1) field multiplications. 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 { 12 // Replace zero Z with one to keep products invertible 13 - var zs: [n]Fe26 = undefined; 14 for (0..n) |i| { 15 - zs[i] = if (points[i].z.isZero()) Fe26.one else points[i].z; 16 } 17 18 // Forward: accumulate Z products 19 - var products: [n]Fe26 = undefined; 20 products[0] = zs[0]; 21 for (1..n) |i| { 22 products[i] = products[i - 1].mul(zs[i]); ··· 27 28 // Backward: recover individual Z inverses and convert to affine 29 // For Jacobian: x = X * zinv², y = Y * zinv³ 30 - var result: [n]AffinePoint26 = undefined; 31 var i: usize = n - 1; 32 while (true) { 33 const z_inv = if (i > 0) inv.mul(products[i - 1]) else inv; 34 if (i > 0) inv = inv.mul(zs[i]); 35 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); ··· 52 } 53 54 /// Build a byte-indexed precomputed table for scalar multiplication. 55 - /// table[i][j] = j * 256^i * base, stored as AffinePoint26. 56 /// 32 subtables × 256 entries = full 256-bit scalar coverage. 57 /// table[i][0] is the identity element (unused in lookups). 58 - pub fn buildByteTable(base: Secp256k1) [32][256]AffinePoint26 { 59 @setEvalBranchQuota(100_000_000); 60 61 - // Phase 1: compute all 32*256 points in Jacobian Fe26 62 var flat: [32 * 256]JacobianPoint = undefined; 63 64 - var cur_base_affine = AffinePoint26.fromStdlib(base.affineCoordinates()); 65 for (0..32) |sub| { 66 flat[sub * 256] = JacobianPoint.identity; 67 flat[sub * 256 + 1] = JacobianPoint.fromAffine(cur_base_affine); ··· 82 const affine_flat = batchToAffine(32 * 256, flat); 83 84 // Reshape to [32][256] 85 - var result: [32][256]AffinePoint26 = undefined; 86 for (0..32) |sub| { 87 for (0..256) |j| { 88 result[sub][j] = affine_flat[sub * 256 + j]; ··· 94 test "batchToAffine matches individual conversion" { 95 const G = Secp256k1.basePoint; 96 const g_affine = G.affineCoordinates(); 97 - const g26 = AffinePoint26.fromStdlib(g_affine); 98 const StdAffineCoordinates = @TypeOf(g_affine); 99 100 const j_id = JacobianPoint.identity;

··· 1 const std = @import("std"); 2 const Secp256k1 = std.crypto.ecc.Secp256k1; 3 + const Fe = @import("field.zig").Fe; 4 const jacobian_mod = @import("jacobian.zig"); 5 const JacobianPoint = jacobian_mod.JacobianPoint; 6 + const AffinePoint = jacobian_mod.AffinePoint; 7 8 /// Batch-convert Jacobian points to affine using Montgomery's trick. 9 /// Uses 1 field inversion + 3*(n-1) field multiplications. 10 /// For Jacobian: x_affine = X * Z^{-2}, y_affine = Y * Z^{-3}. 11 + pub fn batchToAffine(comptime n: usize, points: [n]JacobianPoint) [n]AffinePoint { 12 // Replace zero Z with one to keep products invertible 13 + var zs: [n]Fe = undefined; 14 for (0..n) |i| { 15 + zs[i] = if (points[i].z.isZero()) Fe.one else points[i].z; 16 } 17 18 // Forward: accumulate Z products 19 + var products: [n]Fe = undefined; 20 products[0] = zs[0]; 21 for (1..n) |i| { 22 products[i] = products[i - 1].mul(zs[i]); ··· 27 28 // Backward: recover individual Z inverses and convert to affine 29 // For Jacobian: x = X * zinv², y = Y * zinv³ 30 + var result: [n]AffinePoint = undefined; 31 var i: usize = n - 1; 32 while (true) { 33 const z_inv = if (i > 0) inv.mul(products[i - 1]) else inv; 34 if (i > 0) inv = inv.mul(zs[i]); 35 36 if (points[i].z.isZero()) { 37 + result[i] = AffinePoint.identity; 38 } else { 39 const z_inv2 = z_inv.sq(); 40 const z_inv3 = z_inv2.mul(z_inv); ··· 52 } 53 54 /// Build a byte-indexed precomputed table for scalar multiplication. 55 + /// table[i][j] = j * 256^i * base, stored as AffinePoint. 56 /// 32 subtables × 256 entries = full 256-bit scalar coverage. 57 /// table[i][0] is the identity element (unused in lookups). 58 + pub fn buildByteTable(base: Secp256k1) [32][256]AffinePoint { 59 @setEvalBranchQuota(100_000_000); 60 61 + // Phase 1: compute all 32*256 points in Jacobian Fe 62 var flat: [32 * 256]JacobianPoint = undefined; 63 64 + var cur_base_affine = AffinePoint.fromStdlib(base.affineCoordinates()); 65 for (0..32) |sub| { 66 flat[sub * 256] = JacobianPoint.identity; 67 flat[sub * 256 + 1] = JacobianPoint.fromAffine(cur_base_affine); ··· 82 const affine_flat = batchToAffine(32 * 256, flat); 83 84 // Reshape to [32][256] 85 + var result: [32][256]AffinePoint = undefined; 86 for (0..32) |sub| { 87 for (0..256) |j| { 88 result[sub][j] = affine_flat[sub * 256 + j]; ··· 94 test "batchToAffine matches individual conversion" { 95 const G = Secp256k1.basePoint; 96 const g_affine = G.affineCoordinates(); 97 + const g26 = AffinePoint.fromStdlib(g_affine); 98 const StdAffineCoordinates = @TypeOf(g_affine); 99 100 const j_id = JacobianPoint.identity;

+7 -7

src/endo.zig

··· 1 const std = @import("std"); 2 const Secp256k1 = std.crypto.ecc.Secp256k1; 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; 7 8 /// Cube root of unity in GF(p): beta^3 = 1 mod p. 9 /// phi(x, y) = (beta * x, y) is equivalent to scalar multiplication by lambda, 10 /// but costs only 1 field multiply instead of ~65 doublings. 11 - pub const beta = Fe26.fromInt( 12 55594575648329892869085402983802832744385952214688224221778511981742606582254, 13 ); 14 15 /// Apply the endomorphism to an affine point: (beta * x, y). 16 - pub fn phiAffine(p: AffinePoint26) AffinePoint26 { 17 return .{ .x = p.x.mul(beta), .y = p.y }; 18 } 19 ··· 32 test "phi(G) matches lambda * G" { 33 const G = Secp256k1.basePoint; 34 const g_affine = G.affineCoordinates(); 35 - const g26 = AffinePoint26.fromStdlib(g_affine); 36 37 const phi_G = phiAffine(g26); 38 ··· 52 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); 61 }

··· 1 const std = @import("std"); 2 const Secp256k1 = std.crypto.ecc.Secp256k1; 3 + const Fe = @import("field.zig").Fe; 4 const jacobian_mod = @import("jacobian.zig"); 5 + const AffinePoint = jacobian_mod.AffinePoint; 6 const JacobianPoint = jacobian_mod.JacobianPoint; 7 8 /// Cube root of unity in GF(p): beta^3 = 1 mod p. 9 /// phi(x, y) = (beta * x, y) is equivalent to scalar multiplication by lambda, 10 /// but costs only 1 field multiply instead of ~65 doublings. 11 + pub const beta = Fe.fromInt( 12 55594575648329892869085402983802832744385952214688224221778511981742606582254, 13 ); 14 15 /// Apply the endomorphism to an affine point: (beta * x, y). 16 + pub fn phiAffine(p: AffinePoint) AffinePoint { 17 return .{ .x = p.x.mul(beta), .y = p.y }; 18 } 19 ··· 32 test "phi(G) matches lambda * G" { 33 const G = Secp256k1.basePoint; 34 const g_affine = G.affineCoordinates(); 35 + const g26 = AffinePoint.fromStdlib(g_affine); 36 37 const phi_G = phiAffine(g26); 38 ··· 52 53 // compare via bytes 54 const phi_x = phi_G.x.normalize().toBytes(.big); 55 + const lambda_x = Fe.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 = Fe.fromStdlib(lambda_G_affine.y).normalize().toBytes(.big); 60 try std.testing.expectEqualSlices(u8, &lambda_y, &phi_y); 61 }

+323 -460

src/field.zig

··· 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]; ··· 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 }; ··· 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); ··· 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); ··· 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); ··· 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); ··· 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); ··· 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); ··· 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));

··· 2 const Secp256k1 = std.crypto.ecc.Secp256k1; 3 const StdFe = Secp256k1.Fe; 4 5 + // 5×52-bit unsaturated limb representation (libsecp256k1 style). 6 + // Each limb holds up to 52 bits of value with 12 bits of headroom in u64. 7 + // p = 2^256 - 2^32 - 977, so R = 2^256 mod p = 0x1000003D1. 8 + const M: u64 = 0xFFFFFFFFFFFFF; // 52-bit mask 9 + const R: u64 = 0x1000003D10; // 2^256 mod p (shifted by 4 for limb4's 48-bit width) 10 11 + // secp256k1 prime p in 5×52 12 + const P0: u64 = 0xFFFFEFFFFFC2F; 13 + const P4: u64 = 0xFFFFFFFFFFFF; // 48 bits (limb 4 is only 48 bits wide) 14 15 + /// Optimized secp256k1 field element: 5 × 52-bit limbs in 64-bit words. 16 /// 17 + /// 12 bits of headroom per limb (4 bits in MSW) enables carry-free add/sub 18 + /// for magnitudes up to ~4096. Mul/sq use libsecp256k1's interleaved 19 + /// reduction with R = 2^256 mod p, producing fully reduced output. 20 + pub const Fe = struct { 21 + n: [5]u64, 22 23 + pub const zero: Fe = .{ .n = .{ 0, 0, 0, 0, 0 } }; 24 + pub const one: Fe = .{ .n = .{ 1, 0, 0, 0, 0 } }; 25 26 + /// Pack a 32-byte big-endian value into 5×52-bit limbs. 27 + pub fn fromBytes(b: [32]u8, comptime endian: std.builtin.Endian) Fe { 28 const s = if (endian == .big) b else blk: { 29 var rev: [32]u8 = undefined; 30 for (0..32) |i| rev[i] = b[31 - i]; 31 break :blk rev; 32 }; 33 + // s is big-endian: s[0] is MSB, s[31] is LSB 34 + // limb 0 = bits 0..51 (least significant) 35 + // limb 4 = bits 208..255 (most significant, 48 bits) 36 return .{ .n = .{ 37 + @as(u64, s[31]) | @as(u64, s[30]) << 8 | @as(u64, s[29]) << 16 | 38 + @as(u64, s[28]) << 24 | @as(u64, s[27]) << 32 | @as(u64, s[26]) << 40 | 39 + (@as(u64, s[25]) & 0xF) << 48, 40 + @as(u64, s[25]) >> 4 | @as(u64, s[24]) << 4 | @as(u64, s[23]) << 12 | 41 + @as(u64, s[22]) << 20 | @as(u64, s[21]) << 28 | @as(u64, s[20]) << 36 | 42 + @as(u64, s[19]) << 44, 43 + @as(u64, s[18]) | @as(u64, s[17]) << 8 | @as(u64, s[16]) << 16 | 44 + @as(u64, s[15]) << 24 | @as(u64, s[14]) << 32 | @as(u64, s[13]) << 40 | 45 + (@as(u64, s[12]) & 0xF) << 48, 46 + @as(u64, s[12]) >> 4 | @as(u64, s[11]) << 4 | @as(u64, s[10]) << 12 | 47 + @as(u64, s[9]) << 20 | @as(u64, s[8]) << 28 | @as(u64, s[7]) << 36 | 48 + @as(u64, s[6]) << 44, 49 + @as(u64, s[5]) | @as(u64, s[4]) << 8 | @as(u64, s[3]) << 16 | 50 + @as(u64, s[2]) << 24 | @as(u64, s[1]) << 32 | @as(u64, s[0]) << 40, 51 } }; 52 } 53 54 + /// Unpack 5×52-bit limbs to a 32-byte big-endian value. 55 /// The field value MUST be normalized first. 56 + pub fn toBytes(self: Fe, comptime endian: std.builtin.Endian) [32]u8 { 57 const n = self.n; 58 + var b: [32]u8 = undefined; 59 + // limb 0: bits 0..51 60 b[31] = @truncate(n[0]); 61 b[30] = @truncate(n[0] >> 8); 62 b[29] = @truncate(n[0] >> 16); 63 + b[28] = @truncate(n[0] >> 24); 64 + b[27] = @truncate(n[0] >> 32); 65 + b[26] = @truncate(n[0] >> 40); 66 + // limb 1: bits 52..103 67 + b[25] = @truncate((n[0] >> 48) | (n[1] << 4)); 68 + b[24] = @truncate(n[1] >> 4); 69 + b[23] = @truncate(n[1] >> 12); 70 + b[22] = @truncate(n[1] >> 20); 71 + b[21] = @truncate(n[1] >> 28); 72 + b[20] = @truncate(n[1] >> 36); 73 + b[19] = @truncate(n[1] >> 44); 74 + // limb 2: bits 104..155 75 + b[18] = @truncate(n[2]); 76 + b[17] = @truncate(n[2] >> 8); 77 + b[16] = @truncate(n[2] >> 16); 78 + b[15] = @truncate(n[2] >> 24); 79 + b[14] = @truncate(n[2] >> 32); 80 + b[13] = @truncate(n[2] >> 40); 81 + // limb 3: bits 156..207 82 + b[12] = @truncate((n[2] >> 48) | (n[3] << 4)); 83 + b[11] = @truncate(n[3] >> 4); 84 + b[10] = @truncate(n[3] >> 12); 85 + b[9] = @truncate(n[3] >> 20); 86 + b[8] = @truncate(n[3] >> 28); 87 + b[7] = @truncate(n[3] >> 36); 88 + b[6] = @truncate(n[3] >> 44); 89 + // limb 4: bits 208..255 (48 bits) 90 + b[5] = @truncate(n[4]); 91 + b[4] = @truncate(n[4] >> 8); 92 + b[3] = @truncate(n[4] >> 16); 93 + b[2] = @truncate(n[4] >> 24); 94 + b[1] = @truncate(n[4] >> 32); 95 + b[0] = @truncate(n[4] >> 40); 96 if (endian == .little) { 97 var rev: [32]u8 = undefined; 98 for (0..32) |i| rev[i] = b[31 - i]; ··· 101 return b; 102 } 103 104 + /// Carry propagation + mod-p reduction. 105 + pub fn normalize(self: Fe) Fe { 106 + var t0 = self.n[0]; 107 + var t1 = self.n[1]; 108 + var t2 = self.n[2]; 109 + var t3 = self.n[3]; 110 + var t4 = self.n[4]; 111 112 + // propagate carries 113 + var m = t0 >> 52; 114 + t0 &= M; 115 + t1 += m; 116 + m = t1 >> 52; 117 + t1 &= M; 118 + t2 += m; 119 + m = t2 >> 52; 120 + t2 &= M; 121 + t3 += m; 122 + m = t3 >> 52; 123 + t3 &= M; 124 + t4 += m; 125 + 126 + // reduce: if t4 overflows 48 bits, fold back via p 127 + m = t4 >> 48; 128 + t4 &= 0xFFFFFFFFFFFF; 129 + t0 += m * 0x1000003D1; 130 + // re-propagate 131 + m = t0 >> 52; 132 + t0 &= M; 133 + t1 += m; 134 + m = t1 >> 52; 135 + t1 &= M; 136 + t2 += m; 137 + m = t2 >> 52; 138 + t2 &= M; 139 + t3 += m; 140 + m = t3 >> 52; 141 + t3 &= M; 142 + t4 += m; 143 144 + // conditional subtract p 145 + // check if value >= p 146 + const ge: u64 = if (t4 == 0xFFFFFFFFFFFF and 147 + t3 == M and t2 == M and t1 == M and 148 + t0 >= P0) 1 else 0; 149 + t0 += ge * 0x1000003D1; 150 + m = t0 >> 52; 151 + t0 &= M; 152 + t1 += m; 153 + m = t1 >> 52; 154 + t1 &= M; 155 + t2 += m; 156 + m = t2 >> 52; 157 + t2 &= M; 158 + t3 += m; 159 + m = t3 >> 52; 160 + t3 &= M; 161 + t4 += m; 162 + t4 &= 0xFFFFFFFFFFFF; 163 + 164 + return .{ .n = .{ t0, t1, t2, t3, t4 } }; 165 } 166 167 /// Returns true if the field value equals zero (after normalizing). 168 + pub fn isZero(self: Fe) bool { 169 const f = self.normalize(); 170 + return (f.n[0] | f.n[1] | f.n[2] | f.n[3] | f.n[4]) == 0; 171 } 172 173 /// Returns true if two field values are equal (after normalizing both). 174 + pub fn equivalent(self: Fe, other: Fe) bool { 175 const a = self.normalize(); 176 const b = other.normalize(); 177 + return ((a.n[0] ^ b.n[0]) | (a.n[1] ^ b.n[1]) | (a.n[2] ^ b.n[2]) | 178 + (a.n[3] ^ b.n[3]) | (a.n[4] ^ b.n[4])) == 0; 179 } 180 181 + /// Carry-free addition. 182 + pub fn add(self: Fe, other: Fe) Fe { 183 return .{ .n = .{ 184 + self.n[0] + other.n[0], 185 + self.n[1] + other.n[1], 186 + self.n[2] + other.n[2], 187 + self.n[3] + other.n[3], 188 + self.n[4] + other.n[4], 189 } }; 190 } 191 192 + /// Negate: result = (magnitude+1)*p - self. 193 + pub fn neg(self: Fe, magnitude: u32) Fe { 194 + const m: u64 = @as(u64, magnitude) + 1; 195 return .{ .n = .{ 196 + m * P0 - self.n[0], 197 + m * M - self.n[1], 198 + m * M - self.n[2], 199 + m * M - self.n[3], 200 + m * P4 - self.n[4], 201 } }; 202 } 203 204 + /// sub = self + (-other). 205 + pub fn subMag(self: Fe, other: Fe, other_mag: u32) Fe { 206 return self.add(other.neg(other_mag)); 207 } 208 209 + /// sub with auto-normalize. 210 + pub fn sub(self: Fe, other: Fe) Fe { 211 return self.add(other.normalize().neg(1)); 212 } 213 214 /// Double via carry-free addition. 215 + pub fn dbl(self: Fe) Fe { 216 return self.add(self); 217 } 218 219 + /// Multiply by small integer. 220 + pub fn mulInt(self: Fe, val: u8) Fe { 221 + const v: u64 = val; 222 return .{ .n = .{ 223 + self.n[0] * v, self.n[1] * v, self.n[2] * v, 224 + self.n[3] * v, self.n[4] * v, 225 } }; 226 } 227 228 + /// 5×52 multiplication with interleaved secp256k1 reduction. 229 + /// Ported from libsecp256k1's secp256k1_fe_mul_inner. 230 + /// Uses u128 accumulators with R = 0x1000003D10. 231 + pub fn mul(self: Fe, other: Fe) Fe { 232 const a = self.n; 233 const b = other.n; 234 235 + // column 3: d = p3 236 + var d: u128 = @as(u128, a[0]) * b[3] + @as(u128, a[1]) * b[2] + 237 + @as(u128, a[2]) * b[1] + @as(u128, a[3]) * b[0]; 238 239 + // column 8: c = p8 240 + var c: u128 = @as(u128, a[4]) * b[4]; 241 242 + // reduce p8 into d: d += low64(c) * R, c >>= 64 243 + d += @as(u128, R) * @as(u64, @truncate(c)); 244 + c >>= 64; 245 246 + const t3 = @as(u64, @truncate(d)) & M; 247 + d >>= 52; 248 249 + // column 4: d += p4 250 + d += @as(u128, a[0]) * b[4] + @as(u128, a[1]) * b[3] + 251 + @as(u128, a[2]) * b[2] + @as(u128, a[3]) * b[1] + @as(u128, a[4]) * b[0]; 252 + // reduce remaining c 253 + d += @as(u128, R << 12) * @as(u64, @truncate(c)); 254 255 + var t4 = @as(u64, @truncate(d)) & M; 256 + d >>= 52; 257 + const tx = t4 >> 48; 258 + t4 &= (M >> 4); 259 260 + // column 0: c = p0 261 + c = @as(u128, a[0]) * b[0]; 262 263 + // column 5: d += p5 264 + d += @as(u128, a[1]) * b[4] + @as(u128, a[2]) * b[3] + 265 + @as(u128, a[3]) * b[2] + @as(u128, a[4]) * b[1]; 266 + var u_0 = @as(u64, @truncate(d)) & M; 267 + d >>= 52; 268 269 + u_0 = (u_0 << 4) | tx; 270 + c += @as(u128, u_0) * (R >> 4); 271 272 + const r0 = @as(u64, @truncate(c)) & M; 273 + c >>= 52; 274 275 + // column 1: c += p1 276 + c += @as(u128, a[0]) * b[1] + @as(u128, a[1]) * b[0]; 277 278 + // column 6: d += p6 279 + d += @as(u128, a[2]) * b[4] + @as(u128, a[3]) * b[3] + @as(u128, a[4]) * b[2]; 280 + // reduce low 52 bits of d into c 281 + c += @as(u128, R) * (@as(u64, @truncate(d)) & M); 282 + d >>= 52; 283 284 + const r1 = @as(u64, @truncate(c)) & M; 285 + c >>= 52; 286 287 + // column 2: c += p2 288 + c += @as(u128, a[0]) * b[2] + @as(u128, a[1]) * b[1] + @as(u128, a[2]) * b[0]; 289 290 + // column 7: d += p7 291 + d += @as(u128, a[3]) * b[4] + @as(u128, a[4]) * b[3]; 292 + // reduce: c += low64(d) * R, d >>= 64 293 + c += @as(u128, R) * @as(u64, @truncate(d)); 294 + d >>= 64; 295 296 + const r2 = @as(u64, @truncate(c)) & M; 297 + c >>= 52; 298 299 + // finalize r3, r4 300 + c += @as(u128, R << 12) * @as(u64, @truncate(d)); 301 + c += t3; 302 303 + const r3 = @as(u64, @truncate(c)) & M; 304 + c >>= 52; 305 306 + const r4 = @as(u64, @truncate(c)) + t4; 307 308 + return .{ .n = .{ r0, r1, r2, r3, r4 } }; 309 } 310 311 + /// 5×52 squaring with interleaved secp256k1 reduction. 312 + /// Ported from libsecp256k1's secp256k1_fe_sqr_inner. 313 + pub fn sq(self: Fe) Fe { 314 const a = self.n; 315 316 + // column 3: d = p3 = 2*a0*a3 + 2*a1*a2 317 + var d: u128 = @as(u128, a[0] * 2) * a[3] + @as(u128, a[1] * 2) * a[2]; 318 319 + // column 8: c = p8 = a4*a4 320 + var c: u128 = @as(u128, a[4]) * a[4]; 321 322 + d += @as(u128, R) * @as(u64, @truncate(c)); 323 + c >>= 64; 324 325 + const t3 = @as(u64, @truncate(d)) & M; 326 + d >>= 52; 327 328 + // column 4: d += p4 = a0*2*a4 + 2*a1*a3 + a2*a2 329 + const a4x2 = a[4] * 2; 330 + d += @as(u128, a[0]) * a4x2 + @as(u128, a[1] * 2) * a[3] + @as(u128, a[2]) * a[2]; 331 + d += @as(u128, R << 12) * @as(u64, @truncate(c)); 332 333 + var t4 = @as(u64, @truncate(d)) & M; 334 + d >>= 52; 335 + const tx = t4 >> 48; 336 + t4 &= (M >> 4); 337 338 + // column 0: c = p0 = a0*a0 339 + c = @as(u128, a[0]) * a[0]; 340 341 + // column 5: d += p5 = a1*2*a4 + 2*a2*a3 342 + d += @as(u128, a[1]) * a4x2 + @as(u128, a[2] * 2) * a[3]; 343 + var u_0 = @as(u64, @truncate(d)) & M; 344 + d >>= 52; 345 346 + u_0 = (u_0 << 4) | tx; 347 + c += @as(u128, u_0) * (R >> 4); 348 349 + const r0 = @as(u64, @truncate(c)) & M; 350 + c >>= 52; 351 352 + // column 1: c += p1 = 2*a0*a1 353 + const a0x2 = a[0] * 2; 354 + c += @as(u128, a0x2) * a[1]; 355 356 + // column 6: d += p6 = a2*2*a4 + a3*a3 357 + d += @as(u128, a[2]) * a4x2 + @as(u128, a[3]) * a[3]; 358 + c += @as(u128, R) * (@as(u64, @truncate(d)) & M); 359 + d >>= 52; 360 361 + const r1 = @as(u64, @truncate(c)) & M; 362 + c >>= 52; 363 364 + // column 2: c += p2 = 2*a0*a2 + a1*a1 365 + c += @as(u128, a0x2) * a[2] + @as(u128, a[1]) * a[1]; 366 367 + // column 7: d += p7 = a3*2*a4 368 + d += @as(u128, a[3]) * a4x2; 369 + c += @as(u128, R) * @as(u64, @truncate(d)); 370 + d >>= 64; 371 372 + const r2 = @as(u64, @truncate(c)) & M; 373 + c >>= 52; 374 375 + c += @as(u128, R << 12) * @as(u64, @truncate(d)); 376 + c += t3; 377 378 + const r3 = @as(u64, @truncate(c)) & M; 379 + c >>= 52; 380 381 + const r4 = @as(u64, @truncate(c)) + t4; 382 383 + return .{ .n = .{ r0, r1, r2, r3, r4 } }; 384 } 385 386 + /// Comptime conversion from integer literal to Fe. 387 + pub fn fromInt(comptime val: u256) Fe { 388 var b: [32]u8 = undefined; 389 std.mem.writeInt(u256, &b, val, .big); 390 return fromBytes(b, .big); 391 } 392 393 /// Modular inverse via Fermat's little theorem: a^(p-2) mod p. 394 + pub fn invert(self: Fe) Fe { 395 const a = self; 396 + const a2 = a.sq().mul(a); 397 + const a3 = a2.sq().mul(a); 398 + const a6 = sqn(a3, 3).mul(a3); 399 + const a9 = sqn(a6, 3).mul(a3); 400 + const a11 = sqn(a9, 2).mul(a2); 401 + const a22 = sqn(a11, 11).mul(a11); 402 + const a44 = sqn(a22, 22).mul(a22); 403 + const a88 = sqn(a44, 44).mul(a44); 404 + const a176 = sqn(a88, 88).mul(a88); 405 + const a220 = sqn(a176, 44).mul(a44); 406 + const a223 = sqn(a220, 3).mul(a3); 407 408 + var r = sqn(a223, 23).mul(a22); 409 + r = sqn(r, 5).mul(a); 410 + r = sqn(r, 3).mul(a2); 411 + r = sqn(r, 2).mul(a); 412 return r; 413 } 414 415 + fn sqn(f: Fe, comptime n: usize) Fe { 416 var r = f; 417 for (0..n) |_| r = r.sq(); 418 return r; 419 } 420 421 + /// Convert from stdlib Fe to Fe via byte round-trip. 422 + pub fn fromStdlib(fe: StdFe) Fe { 423 return fromBytes(fe.toBytes(.big), .big); 424 } 425 426 /// Convert to stdlib Fe via byte round-trip. 427 + pub fn toStdlib(self: Fe) StdFe { 428 return StdFe.fromBytes(self.normalize().toBytes(.big), .big) catch unreachable; 429 } 430 }; ··· 444 } 445 446 test "fromBytes/toBytes round-trip" { 447 const gx_bytes = Secp256k1.basePoint.affineCoordinates().x.toBytes(.big); 448 + const fe = Fe.fromBytes(gx_bytes, .big); 449 const out = fe.normalize().toBytes(.big); 450 try testing.expectEqualSlices(u8, &gx_bytes, &out); 451 } 452 453 test "fromBytes/toBytes little-endian round-trip" { 454 const gx_bytes = Secp256k1.basePoint.affineCoordinates().x.toBytes(.little); 455 + const fe = Fe.fromBytes(gx_bytes, .little); 456 const out = fe.normalize().toBytes(.little); 457 try testing.expectEqualSlices(u8, &gx_bytes, &out); 458 } 459 460 test "zero and one" { 461 + try testing.expect(Fe.zero.isZero()); 462 + try testing.expect(!Fe.one.isZero()); 463 + try testing.expect(Fe.one.equivalent(Fe.one)); 464 + try testing.expect(!Fe.zero.equivalent(Fe.one)); 465 } 466 467 test "add matches stdlib" { 468 const g = Secp256k1.basePoint.affineCoordinates(); 469 + const a = Fe.fromStdlib(g.x); 470 + const b = Fe.fromStdlib(g.y); 471 const result = a.add(b).normalize().toBytes(.big); 472 const expected = g.x.add(g.y).toBytes(.big); 473 try testing.expectEqualSlices(u8, &expected, &result); ··· 475 476 test "sub matches stdlib" { 477 const g = Secp256k1.basePoint.affineCoordinates(); 478 + const a = Fe.fromStdlib(g.x); 479 + const b = Fe.fromStdlib(g.y); 480 const result = a.sub(b).normalize().toBytes(.big); 481 const expected = g.x.sub(g.y).toBytes(.big); 482 try testing.expectEqualSlices(u8, &expected, &result); ··· 484 485 test "neg matches stdlib" { 486 const g = Secp256k1.basePoint.affineCoordinates(); 487 + const a = Fe.fromStdlib(g.x); 488 const result = a.neg(1).normalize().toBytes(.big); 489 const expected = StdFe.zero.sub(g.x).toBytes(.big); 490 try testing.expectEqualSlices(u8, &expected, &result); ··· 492 493 test "mul matches stdlib" { 494 const g = Secp256k1.basePoint.affineCoordinates(); 495 + const a = Fe.fromStdlib(g.x); 496 + const b = Fe.fromStdlib(g.y); 497 const result = a.mul(b).normalize().toBytes(.big); 498 const expected = stdFeMul(g.x, g.y).toBytes(.big); 499 try testing.expectEqualSlices(u8, &expected, &result); ··· 501 502 test "sq matches stdlib" { 503 const g = Secp256k1.basePoint.affineCoordinates(); 504 + const a = Fe.fromStdlib(g.x); 505 const result = a.sq().normalize().toBytes(.big); 506 const expected = stdFeSq(g.x).toBytes(.big); 507 try testing.expectEqualSlices(u8, &expected, &result); ··· 509 510 test "sq equals mul(a, a)" { 511 const g = Secp256k1.basePoint.affineCoordinates(); 512 + const a = Fe.fromStdlib(g.x); 513 try testing.expect(a.sq().equivalent(a.mul(a))); 514 515 + const b = Fe.fromStdlib(g.y); 516 try testing.expect(b.sq().equivalent(b.mul(b))); 517 } 518 519 test "mul(a, invert(a)) == one" { 520 const g = Secp256k1.basePoint.affineCoordinates(); 521 + const a = Fe.fromStdlib(g.x); 522 const inv = a.invert(); 523 const product = a.mul(inv); 524 + try testing.expect(product.equivalent(Fe.one)); 525 } 526 527 test "mul(a, invert(a)) == one for y" { 528 const g = Secp256k1.basePoint.affineCoordinates(); 529 + const a = Fe.fromStdlib(g.y); 530 const inv = a.invert(); 531 const product = a.mul(inv); 532 + try testing.expect(product.equivalent(Fe.one)); 533 } 534 535 test "edge: p-1" { 536 + const p_minus_1 = Fe.fromInt( 537 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2E, 538 ); 539 try testing.expect(!p_minus_1.isZero()); 540 + try testing.expect(p_minus_1.add(Fe.one).isZero()); 541 } 542 543 test "fromInt constants" { 544 + const beta = Fe.fromInt( 545 55594575648329892869085402983802832744385952214688224221778511981742606582254, 546 ); 547 try testing.expect(!beta.isZero()); 548 const beta_cubed = beta.sq().mul(beta); 549 + try testing.expect(beta_cubed.equivalent(Fe.one)); 550 } 551 552 test "dbl matches add(self)" { 553 const g = Secp256k1.basePoint.affineCoordinates(); 554 + const a = Fe.fromStdlib(g.x); 555 try testing.expect(a.dbl().equivalent(a.add(a))); 556 } 557 558 test "mulInt matches repeated add" { 559 const g = Secp256k1.basePoint.affineCoordinates(); 560 + const a = Fe.fromStdlib(g.x); 561 try testing.expect(a.mulInt(3).equivalent(a.add(a).add(a))); 562 } 563 564 test "multiple mul/sq operations stay consistent" { 565 const g = Secp256k1.basePoint.affineCoordinates(); 566 + const gx26 = Fe.fromStdlib(g.x); 567 + const gy26 = Fe.fromStdlib(g.y); 568 569 const result = gx26.sq().mul(gy26); 570 const expected = stdFeSq(g.x).mul(g.y); ··· 573 574 test "normalize idempotent" { 575 const g = Secp256k1.basePoint.affineCoordinates(); 576 + const a = Fe.fromStdlib(g.x); 577 const n1 = a.normalize(); 578 const n2 = n1.normalize(); 579 try testing.expectEqualSlices(u8, &std.mem.toBytes(n1.n), &std.mem.toBytes(n2.n));

+31 -31

src/jacobian.zig

··· 1 const std = @import("std"); 2 const Secp256k1 = std.crypto.ecc.Secp256k1; 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 }; 13 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³). 33 /// Uses a=0 specialized formulas for fewer field operations than 34 /// the stdlib's complete projective formulas. 35 pub const JacobianPoint = struct { 36 - x: Fe26, 37 - y: Fe26, 38 - z: Fe26, 39 40 - pub const identity: JacobianPoint = .{ .x = Fe26.zero, .y = Fe26.one, .z = Fe26.zero }; 41 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 }; 45 } 46 47 - /// Convert stdlib affine to Jacobian via Fe26. 48 pub fn fromStdlibAffine(ac: StdAffineCoordinates) JacobianPoint { 49 - return fromAffine(AffinePoint26.fromStdlib(ac)); 50 } 51 52 /// Doubling for a=0 curves (secp256k1). 2M + 5S + 3 normalize. ··· 71 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 { 75 if (self.z.isZero()) return fromAffine(q); 76 77 const z1z1 = self.z.sq(); // mag 1 ··· 98 } 99 100 /// Mixed subtraction: Jacobian - Affine → Jacobian. 101 - pub fn subMixed(self: JacobianPoint, q: AffinePoint26) JacobianPoint { 102 return self.addMixed(q.neg()); 103 } 104 ··· 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(); ··· 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); ··· 182 }; 183 184 // ============================================================ 185 - // tests — verify Fe26 Jacobian matches stdlib results 186 // ============================================================ 187 188 fn toStdlibAffine(j: JacobianPoint) StdAffineCoordinates { ··· 205 test "addMixed matches stdlib add" { 206 const G = Secp256k1.basePoint; 207 const g_affine = G.affineCoordinates(); 208 - const g26 = AffinePoint26.fromStdlib(g_affine); 209 210 // 2G + G = 3G 211 const j_3g = JacobianPoint.fromAffine(g26).dbl().addMixed(g26); ··· 220 test "identity + affine = affine" { 221 const G = Secp256k1.basePoint; 222 const g_affine = G.affineCoordinates(); 223 - const g26 = AffinePoint26.fromStdlib(g_affine); 224 225 const result = JacobianPoint.identity.addMixed(g26); 226 const result_affine = toStdlibAffine(result); ··· 232 test "subMixed is inverse of addMixed" { 233 const G = Secp256k1.basePoint; 234 const g_affine = G.affineCoordinates(); 235 - const g26 = AffinePoint26.fromStdlib(g_affine); 236 237 // 3G - G = 2G 238 const j_3g = JacobianPoint.fromAffine(g26).dbl().addMixed(g26); ··· 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(); ··· 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);

··· 1 const std = @import("std"); 2 const Secp256k1 = std.crypto.ecc.Secp256k1; 3 + const Fe = @import("field.zig").Fe; 4 const StdFe = Secp256k1.Fe; 5 const StdAffineCoordinates = @TypeOf(Secp256k1.basePoint.affineCoordinates()); 6 7 + /// Affine point in Fe representation. 8 + pub const AffinePoint = struct { 9 + x: Fe, 10 + y: Fe, 11 12 + pub const identity: AffinePoint = .{ .x = Fe.zero, .y = Fe.zero }; 13 14 + pub fn neg(self: AffinePoint) AffinePoint { 15 return .{ .x = self.x, .y = self.y.neg(1) }; 16 } 17 18 + pub fn fromStdlib(ac: StdAffineCoordinates) AffinePoint { 19 return .{ 20 + .x = Fe.fromBytes(ac.x.toBytes(.big), .big), 21 + .y = Fe.fromBytes(ac.y.toBytes(.big), .big), 22 }; 23 } 24 25 + pub fn fromStdlibIdentity(ac: StdAffineCoordinates) AffinePoint { 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 Fe: affine (X/Z², Y/Z³). 33 /// Uses a=0 specialized formulas for fewer field operations than 34 /// the stdlib's complete projective formulas. 35 pub const JacobianPoint = struct { 36 + x: Fe, 37 + y: Fe, 38 + z: Fe, 39 40 + pub const identity: JacobianPoint = .{ .x = Fe.zero, .y = Fe.one, .z = Fe.zero }; 41 42 + /// Convert Fe affine (x, y) to Jacobian (x, y, 1). 43 + pub fn fromAffine(p: AffinePoint) JacobianPoint { 44 + return .{ .x = p.x, .y = p.y, .z = Fe.one }; 45 } 46 47 + /// Convert stdlib affine to Jacobian via Fe. 48 pub fn fromStdlibAffine(ac: StdAffineCoordinates) JacobianPoint { 49 + return fromAffine(AffinePoint.fromStdlib(ac)); 50 } 51 52 /// Doubling for a=0 curves (secp256k1). 2M + 5S + 3 normalize. ··· 71 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: AffinePoint) JacobianPoint { 75 if (self.z.isZero()) return fromAffine(q); 76 77 const z1z1 = self.z.sq(); // mag 1 ··· 98 } 99 100 /// Mixed subtraction: Jacobian - Affine → Jacobian. 101 + pub fn subMixed(self: JacobianPoint, q: AffinePoint) JacobianPoint { 102 return self.addMixed(q.neg()); 103 } 104 ··· 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: Fe) bool { 150 // Jacobian: x_affine = X / Z² 151 // So check: r * Z² == X 152 const z2 = self.z.sq(); ··· 155 156 // rare overflow case: x_affine could be in [n, p) 157 // check (r + n) * Z² == X 158 + const n_fe26 = comptime Fe.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); ··· 182 }; 183 184 // ============================================================ 185 + // tests — verify Fe Jacobian matches stdlib results 186 // ============================================================ 187 188 fn toStdlibAffine(j: JacobianPoint) StdAffineCoordinates { ··· 205 test "addMixed matches stdlib add" { 206 const G = Secp256k1.basePoint; 207 const g_affine = G.affineCoordinates(); 208 + const g26 = AffinePoint.fromStdlib(g_affine); 209 210 // 2G + G = 3G 211 const j_3g = JacobianPoint.fromAffine(g26).dbl().addMixed(g26); ··· 220 test "identity + affine = affine" { 221 const G = Secp256k1.basePoint; 222 const g_affine = G.affineCoordinates(); 223 + const g26 = AffinePoint.fromStdlib(g_affine); 224 225 const result = JacobianPoint.identity.addMixed(g26); 226 const result_affine = toStdlibAffine(result); ··· 232 test "subMixed is inverse of addMixed" { 233 const G = Secp256k1.basePoint; 234 const g_affine = G.affineCoordinates(); 235 + const g26 = AffinePoint.fromStdlib(g_affine); 236 237 // 3G - G = 2G 238 const j_3g = JacobianPoint.fromAffine(g26).dbl().addMixed(g26); ··· 262 test "full Jacobian add matches stdlib" { 263 const G = Secp256k1.basePoint; 264 const g_affine = G.affineCoordinates(); 265 + const g26 = AffinePoint.fromStdlib(g_affine); 266 267 // 2G + 3G = 5G 268 const j_2g = JacobianPoint.fromAffine(g26).dbl(); ··· 280 test "add with identity" { 281 const G = Secp256k1.basePoint; 282 const g_affine = G.affineCoordinates(); 283 + const g26 = AffinePoint.fromStdlib(g_affine); 284 285 const j_g = JacobianPoint.fromAffine(g26); 286 const result = j_g.add(JacobianPoint.identity);

+5 -5

src/point.zig

··· 1 const std = @import("std"); 2 const jacobian_mod = @import("jacobian.zig"); 3 const JacobianPoint = jacobian_mod.JacobianPoint; 4 - const AffinePoint26 = jacobian_mod.AffinePoint26; 5 const endo = @import("endo.zig"); 6 const affine_mod = @import("affine.zig"); 7 const Secp256k1 = std.crypto.ecc.Secp256k1; 8 9 /// Build a precomputation table: [identity, p, 2p, 3p, ..., count*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); ··· 21 22 /// Apply the endomorphism to each entry in an affine precompute table. 23 /// phi(n*P) = n*phi(P), so transforming the table gives a valid table for phi(P). 24 - pub fn phiTableAffine(comptime n: usize, table: [n]AffinePoint26) [n]AffinePoint26 { 25 - var result: [n]AffinePoint26 = undefined; 26 for (0..n) |i| { 27 result[i] = endo.phiAffine(table[i]); 28 } ··· 60 test "precompute table correctness" { 61 const G = Secp256k1.basePoint; 62 const g_affine = G.affineCoordinates(); 63 - const g26 = AffinePoint26.fromStdlib(g_affine); 64 65 const table = precompute(g26, 8); 66

··· 1 const std = @import("std"); 2 const jacobian_mod = @import("jacobian.zig"); 3 const JacobianPoint = jacobian_mod.JacobianPoint; 4 + const AffinePoint = jacobian_mod.AffinePoint; 5 const endo = @import("endo.zig"); 6 const affine_mod = @import("affine.zig"); 7 const Secp256k1 = std.crypto.ecc.Secp256k1; 8 9 /// Build a precomputation table: [identity, p, 2p, 3p, ..., count*p]. 10 /// Uses Jacobian arithmetic for efficient computation. 11 + pub fn precompute(p: AffinePoint, comptime count: usize) [1 + count]JacobianPoint { 12 var pc: [1 + count]JacobianPoint = undefined; 13 pc[0] = JacobianPoint.identity; 14 pc[1] = JacobianPoint.fromAffine(p); ··· 21 22 /// Apply the endomorphism to each entry in an affine precompute table. 23 /// phi(n*P) = n*phi(P), so transforming the table gives a valid table for phi(P). 24 + pub fn phiTableAffine(comptime n: usize, table: [n]AffinePoint) [n]AffinePoint { 25 + var result: [n]AffinePoint = undefined; 26 for (0..n) |i| { 27 result[i] = endo.phiAffine(table[i]); 28 } ··· 60 test "precompute table correctness" { 61 const G = Secp256k1.basePoint; 62 const g_affine = G.affineCoordinates(); 63 + const g26 = AffinePoint.fromStdlib(g_affine); 64 65 const table = precompute(g26, 8); 66

+329 -21

src/verify.zig

··· 1 const std = @import("std"); 2 const Secp256k1 = std.crypto.ecc.Secp256k1; 3 - const Fe26 = @import("field.zig").Fe26; 4 const scalar = Secp256k1.scalar; 5 6 const endo = @import("endo.zig"); ··· 8 const affine_mod = @import("affine.zig"); 9 const jacobian_mod = @import("jacobian.zig"); 10 const JacobianPoint = jacobian_mod.JacobianPoint; 11 - const AffinePoint26 = jacobian_mod.AffinePoint26; 12 13 pub const VerifyError = std.crypto.errors.IdentityElementError || 14 std.crypto.errors.NonCanonicalError || 15 error{SignatureVerificationFailed}; 16 17 - /// Precomputed base point table: G_TABLE[i][j] = j * 256^i * G in Fe26 affine. 18 /// 32 subtables × 256 entries = 8192 affine points. 19 /// Enables u1*G computation via byte-at-a-time lookup with zero doublings. 20 - const G_TABLE: [32][256]AffinePoint26 = blk: { 21 - @setEvalBranchQuota(100_000_000); 22 - break :blk affine_mod.buildByteTable(Secp256k1.basePoint); 23 - }; 24 25 /// Reduce a 32-byte big-endian value to a scalar (mod n). 26 fn reduceToScalar(h: [32]u8) scalar.Scalar { ··· 35 /// 1. u1*G via precomputed byte table: ~32 mixed adds, zero doublings 36 /// 2. u2*Q via projective-table GLV: no field inversion, 4-bit windowed 37 /// 3. Jacobian comparison (no field inversion for final check) 38 - /// 4. All arithmetic in Fe26 (10×26-bit) — no Montgomery form overhead 39 pub fn verify(sig_r: [32]u8, sig_s: [32]u8, msg_hash: [32]u8, public_key: Secp256k1) VerifyError!void { 40 // parse and validate r, s 41 const r_sc = scalar.Scalar.fromBytes(sig_r, .big) catch return error.SignatureVerificationFailed; ··· 44 45 // scalar_u1 = z * s^-1, scalar_u2 = r * s^-1 46 const z = reduceToScalar(msg_hash); 47 - const s_inv = s_sc.invert(); 48 const scalar_u1 = z.mul(s_inv).toBytes(.big); 49 const scalar_u2 = r_sc.mul(s_inv).toBytes(.little); 50 ··· 64 split_u2.r2 = scalar.neg(split_u2.r2, .little) catch zero_s; 65 neg_p_phi = true; 66 } 67 - const pk_affine26 = AffinePoint26.fromStdlib(public_key.affineCoordinates()); 68 const r2 = publicKeyMulProjective(split_u2, neg_p, neg_p_phi, pk_affine26); 69 70 // 3. combine results in Jacobian, compare without inversion ··· 74 if (q.z.isZero()) return error.SignatureVerificationFailed; 75 76 // Jacobian comparison: r * Z² == X (mod p) 77 - const r_fe = Fe26.fromBytes(sig_r, .big); 78 if (!q.jacobianCompare(r_fe)) { 79 return error.SignatureVerificationFailed; 80 } ··· 84 /// Direct 32-byte lookup: G_TABLE[i][scalar[i]], no GLV decomposition. 85 /// Cost: ~32 mixed Jacobian-affine additions, zero doublings. 86 fn basePointMul(scalar_u1: [32]u8) JacobianPoint { 87 var acc = JacobianPoint.identity; 88 89 - // G_TABLE[i] corresponds to byte position i (little-endian). 90 // scalar_u1 is big-endian, so byte 31 is least significant. 91 - // G_TABLE[i] ↔ scalar_u1[31-i] 92 for (0..32) |i| { 93 const byte = scalar_u1[31 - i]; 94 if (byte != 0) { 95 - acc = acc.addMixed(G_TABLE[i][byte]); 96 } 97 } 98 ··· 106 split: endo.SplitScalar, 107 neg_p: bool, 108 neg_p_phi: bool, 109 - pk_affine: AffinePoint26, 110 ) JacobianPoint { 111 // Build 9-entry Jacobian tables: [identity, 1P, 2P, ..., 8P] 112 const pk_table = point.precompute(pk_affine, 8); ··· 139 return q; 140 } 141 142 - test "G_TABLE[0][1] is the base point" { 143 const G = Secp256k1.basePoint; 144 const g_affine = G.affineCoordinates(); 145 146 - // G_TABLE[0][1] should be 1*G 147 - const table_g = G_TABLE[0][1]; 148 const table_g_x = table_g.x.toStdlib(); 149 const table_g_y = table_g.y.toStdlib(); 150 ··· 152 try std.testing.expect(table_g_y.equivalent(g_affine.y)); 153 } 154 155 - test "G_TABLE[0][2] is 2G" { 156 const G = Secp256k1.basePoint; 157 const g2 = G.dbl().affineCoordinates(); 158 159 - const table_2g = G_TABLE[0][2]; 160 const table_2g_x = table_2g.x.toStdlib(); 161 const table_2g_y = table_2g.y.toStdlib(); 162 ··· 166 167 test "comptime addMixed gives correct 2G" { 168 @setEvalBranchQuota(100_000_000); 169 - const g_affine = comptime AffinePoint26.fromStdlib(Secp256k1.basePoint.affineCoordinates()); 170 const j_g = comptime JacobianPoint.fromAffine(g_affine); 171 const j_2g = comptime j_g.addMixed(g_affine); 172

··· 1 const std = @import("std"); 2 const Secp256k1 = std.crypto.ecc.Secp256k1; 3 + const Fe = @import("field.zig").Fe; 4 const scalar = Secp256k1.scalar; 5 6 const endo = @import("endo.zig"); ··· 8 const affine_mod = @import("affine.zig"); 9 const jacobian_mod = @import("jacobian.zig"); 10 const JacobianPoint = jacobian_mod.JacobianPoint; 11 + const AffinePoint = jacobian_mod.AffinePoint; 12 13 pub 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). 22 + var g_table_storage: [32][256]AffinePoint = undefined; 23 + var g_table_ready: bool = false; 24 + 25 + fn 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). 34 fn reduceToScalar(h: [32]u8) scalar.Scalar { ··· 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 47 pub 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; ··· 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 ··· 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 ··· 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 } ··· 92 /// Direct 32-byte lookup: G_TABLE[i][scalar[i]], no GLV decomposition. 93 /// Cost: ~32 mixed Jacobian-affine additions, zero doublings. 94 fn 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 ··· 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); ··· 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. 157 + fn 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) 395 + inline 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 + 403 + test "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 + 449 + test "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 ··· 459 try std.testing.expect(table_g_y.equivalent(g_affine.y)); 460 } 461 462 + test "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 ··· 474 475 test "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