diff --git a/quartz/src/main/c/secp256k1/field.c b/quartz/src/main/c/secp256k1/field.c index 6a97dba129..41146764e6 100644 --- a/quartz/src/main/c/secp256k1/field.c +++ b/quartz/src/main/c/secp256k1/field.c @@ -1,220 +1,182 @@ /* * Copyright (c) 2025 Vitor Pamplona * - * Field arithmetic mod p = 2^256 - 2^32 - 977, using 5x52-bit limbs. + * Field arithmetic modulo p = 2^256 - 2^32 - 977 using 4x64-bit limbs. * - * The 5x52 representation gives 12 bits of headroom per limb, enabling - * lazy reduction: multiple adds/subs can chain without normalizing. This - * is the key advantage over the Kotlin 4x64-bit approach which must - * reduce after every single add/sub. - * - * On ARM64/x86_64 with __int128: each limb multiply is a single MUL+UMULH - * (ARM64) or MULQ (x86_64) instruction pair, vs the Kotlin version which - * needs Math.multiplyHigh() + unsigned correction (5 JVM instructions). + * Same representation as the Kotlin Fe4 class. Performance analysis showed + * 4x64 is faster than 5x52 because fewer multiplies (16 vs 25) outweighs + * the lazy reduction advantage of 5x52 on both JVM and native. */ #include "field.h" #include -/* ==================== Field Multiplication ==================== */ +#define FIELD_C 0x1000003D1ULL + +/* ==================== 8-limb product computation ==================== */ /* - * Multiply: r = a * b mod p + * Compute 512-bit product of two 256-bit numbers in 4x64 representation. + * Output: 8 limbs in little-endian. Uses __int128 for 64x64->128 products. * - * Uses schoolbook 5x5 multiplication with __int128 for the 64x64->128 products. - * After computing the full 10-limb product, reduces mod p using: - * 2^260 = 2^4 * (2^32 + 977) = 16 * 0x1000003D1 mod p - * - * The reduction folds the high limbs back using the secp256k1 constant R = 2^256 mod p. - * For 5x52 limbs, the folding constant for limb 5 is 0x1000003D1 (since 2^260 = 2^4 * (2^32+977)). + * The key insight: we CAN'T accumulate multiple 128-bit products in a single + * uint128 because 4 products overflow (4 * 2^128 > 2^128). Instead, we compute + * each column separately and propagate carries explicitly. */ #if HAVE_INT128 -void fe_mul(secp256k1_fe *r, const secp256k1_fe *a, const secp256k1_fe *b) { - const uint64_t M = FE_LIMB_MASK; - /* 2^260 mod p: each limb is 52 bits, so position 5 is at 260 bits. - * 2^260 mod p = 2^4 * (2^256 mod p) = 16 * 0x1000003D1 = 0x10000003D10 */ - const uint64_t R = 0x1000003D10ULL; /* 2^260 mod p = 16 * (2^32 + 977) */ - uint128_t c, d; - uint64_t t0, t1, t2, t3, t4; - uint64_t a0 = a->d[0], a1 = a->d[1], a2 = a->d[2], a3 = a->d[3], a4 = a->d[4]; - uint64_t b0 = b->d[0], b1 = b->d[1], b2 = b->d[2], b3 = b->d[3], b4 = b->d[4]; - +static void mul_wide(uint64_t out[8], const uint64_t a[4], const uint64_t b[4]) { /* - * libsecp256k1-style split R-folding. The sum of folded products c can be up to - * ~106 bits. c*R would overflow uint128 (up to 140 bits). Instead, we split: - * d += (uint64_t)c * R (low 64 bits × R, fits in ~98 bits) - * carry (c >> 64) * R into the next limb's accumulator + * Schoolbook 4x4 multiplication into 8 limbs. Uses a row-based approach: + * multiply each a[i] by the full b[0..3] vector and accumulate into out. + * This avoids the column-based carry overflow problem. */ + uint128_t acc; + uint64_t carry; - /* - * Correct approach: expand c*R into individual products a[i]*b[j]*R. - * Each a[i]*b[j]*R < 2^52 * 2^52 * 2^34 = 2^138 which overflows uint128. - * - * Real solution: split R into a[i]*R (uint128) before multiplying by b[j]. - * a[i]*R is at most 2^86 which fits in uint128. Then (uint128)(a[i]*R) * b[j] - * is at most 2^138 which ALSO overflows uint128! - * - * The ACTUAL libsecp256k1 trick: use d and c as two separate accumulators. - * d accumulates the direct products. c accumulates the folded products with - * a DIFFERENT carry chain. Let me just do the schoolbook 10-limb product - * and then reduce. - */ - { - /* Full 10-limb product, then reduce mod p using R = 2^260 / 2^4 fold */ - uint128_t p[10] = {0}; - int i, j; - for (i = 0; i < 5; i++) { - for (j = 0; j < 5; j++) { - p[i+j] += (uint128_t)a->d[i] * b->d[j]; + /* Row 0: out += a[0] * b */ + acc = (uint128_t)a[0] * b[0]; + out[0] = (uint64_t)acc; acc >>= 64; + acc += (uint128_t)a[0] * b[1]; + out[1] = (uint64_t)acc; acc >>= 64; + acc += (uint128_t)a[0] * b[2]; + out[2] = (uint64_t)acc; acc >>= 64; + acc += (uint128_t)a[0] * b[3]; + out[3] = (uint64_t)acc; + out[4] = (uint64_t)(acc >> 64); + out[5] = out[6] = out[7] = 0; + + /* Row 1: out[1..5] += a[1] * b */ + acc = (uint128_t)out[1] + (uint128_t)a[1] * b[0]; + out[1] = (uint64_t)acc; acc >>= 64; + acc += (uint128_t)out[2] + (uint128_t)a[1] * b[1]; + out[2] = (uint64_t)acc; acc >>= 64; + acc += (uint128_t)out[3] + (uint128_t)a[1] * b[2]; + out[3] = (uint64_t)acc; acc >>= 64; + acc += (uint128_t)out[4] + (uint128_t)a[1] * b[3]; + out[4] = (uint64_t)acc; + out[5] = (uint64_t)(acc >> 64); + + /* Row 2: out[2..6] += a[2] * b */ + acc = (uint128_t)out[2] + (uint128_t)a[2] * b[0]; + out[2] = (uint64_t)acc; acc >>= 64; + acc += (uint128_t)out[3] + (uint128_t)a[2] * b[1]; + out[3] = (uint64_t)acc; acc >>= 64; + acc += (uint128_t)out[4] + (uint128_t)a[2] * b[2]; + out[4] = (uint64_t)acc; acc >>= 64; + acc += (uint128_t)out[5] + (uint128_t)a[2] * b[3]; + out[5] = (uint64_t)acc; + out[6] = (uint64_t)(acc >> 64); + + /* Row 3: out[3..7] += a[3] * b */ + acc = (uint128_t)out[3] + (uint128_t)a[3] * b[0]; + out[3] = (uint64_t)acc; acc >>= 64; + acc += (uint128_t)out[4] + (uint128_t)a[3] * b[1]; + out[4] = (uint64_t)acc; acc >>= 64; + acc += (uint128_t)out[5] + (uint128_t)a[3] * b[2]; + out[5] = (uint64_t)acc; acc >>= 64; + acc += (uint128_t)out[6] + (uint128_t)a[3] * b[3]; + out[6] = (uint64_t)acc; + out[7] = (uint64_t)(acc >> 64); +} + +/* + * Reduce a 512-bit value (8 limbs) modulo p. + * Uses: 2^256 ≡ C (mod p) where C = 0x1000003D1. + * Two reduction rounds: first folds hi[0..3] into lo[0..3] using C, + * second handles any remaining overflow. + */ +static void reduce_wide(secp256k1_fe *r, const uint64_t w[8]) { + uint128_t acc; + + /* Round 1: result = w[0..3] + w[4..7] * C */ + acc = (uint128_t)w[0] + (uint128_t)w[4] * FIELD_C; + r->d[0] = (uint64_t)acc; + acc >>= 64; + + acc += (uint128_t)w[1] + (uint128_t)w[5] * FIELD_C; + r->d[1] = (uint64_t)acc; + acc >>= 64; + + acc += (uint128_t)w[2] + (uint128_t)w[6] * FIELD_C; + r->d[2] = (uint64_t)acc; + acc >>= 64; + + acc += (uint128_t)w[3] + (uint128_t)w[7] * FIELD_C; + r->d[3] = (uint64_t)acc; + uint64_t carry = (uint64_t)(acc >> 64); + + /* Round 2: fold remaining carry */ + if (carry) { + acc = (uint128_t)r->d[0] + (uint128_t)carry * FIELD_C; + r->d[0] = (uint64_t)acc; + carry = (uint64_t)(acc >> 64); + if (carry) { + r->d[1] += carry; + if (r->d[1] < carry) { + r->d[2]++; + if (r->d[2] == 0) r->d[3]++; } } - /* Propagate carries in the 10-limb product */ - for (i = 0; i < 9; i++) { - p[i+1] += p[i] >> 52; - p[i] &= M; - } - /* Fold high limbs using R = 0x1000003D1 (2^256 mod p in 5x52) */ - /* Limbs 5-9 fold into 0-4: p[i+5] * R adds to p[i] */ - for (i = 4; i >= 0; i--) { - if (i + 5 <= 9 && p[i+5]) { - uint128_t fold = p[i+5] * R; - p[i] += fold; - } - } - /* Propagate carries again */ - for (i = 0; i < 4; i++) { - p[i+1] += p[i] >> 52; - p[i] &= M; - } - /* Fold any remaining overflow from limb 4. - * Limb 4 overflow is at bit position 4*52+48 = 256, so use 2^256 mod p = 0x1000003D1 */ - if (p[4] >> 48) { - uint64_t overflow = (uint64_t)(p[4] >> 48); - p[4] &= 0xFFFFFFFFFFFFULL; - p[0] += (uint128_t)overflow * 0x1000003D1ULL; - p[1] += p[0] >> 52; p[0] &= M; - p[2] += p[1] >> 52; p[1] &= M; - p[3] += p[2] >> 52; p[2] &= M; - p[4] += p[3] >> 52; p[3] &= M; - } - r->d[0] = (uint64_t)p[0]; r->d[1] = (uint64_t)p[1]; r->d[2] = (uint64_t)p[2]; - r->d[3] = (uint64_t)p[3]; r->d[4] = (uint64_t)p[4]; } + + fe_normalize(r); +} + +void fe_mul(secp256k1_fe *r, const secp256k1_fe *a, const secp256k1_fe *b) { + uint64_t w[8]; + mul_wide(w, a->d, b->d); + reduce_wide(r, w); } void fe_sqr(secp256k1_fe *r, const secp256k1_fe *a) { - const uint64_t M = FE_LIMB_MASK; - /* 2^260 mod p: each limb is 52 bits, so position 5 is at 260 bits. - * 2^260 mod p = 2^4 * (2^256 mod p) = 16 * 0x1000003D1 = 0x10000003D10 */ - const uint64_t R = 0x1000003D10ULL; /* 2^260 mod p = 16 * (2^32 + 977) */ - uint128_t c, d; - uint64_t t0, t1, t2, t3, t4; - uint64_t a0 = a->d[0], a1 = a->d[1], a2 = a->d[2], a3 = a->d[3], a4 = a->d[4]; - - /* Same split R-folding approach as fe_mul, with doubled cross-products */ - - /* Use schoolbook product + reduction (same as fe_mul but with doubled cross-products) */ - { - uint128_t p[10] = {0}; - int i, j; - for (i = 0; i < 5; i++) { - for (j = 0; j < 5; j++) { - p[i+j] += (uint128_t)a->d[i] * a->d[j]; - } - } - for (i = 0; i < 9; i++) { - p[i+1] += p[i] >> 52; - p[i] &= M; - } - for (i = 4; i >= 0; i--) { - if (i + 5 <= 9 && p[i+5]) { - p[i] += p[i+5] * R; - } - } - for (i = 0; i < 4; i++) { - p[i+1] += p[i] >> 52; - p[i] &= M; - } - if (p[4] >> 48) { - uint64_t overflow = (uint64_t)(p[4] >> 48); - p[4] &= 0xFFFFFFFFFFFFULL; - p[0] += (uint128_t)overflow * 0x1000003D1ULL; /* 2^256 mod p */ - p[1] += p[0] >> 52; p[0] &= M; - p[2] += p[1] >> 52; p[1] &= M; - p[3] += p[2] >> 52; p[2] &= M; - p[4] += p[3] >> 52; p[3] &= M; - } - t0 = (uint64_t)p[0]; t1 = (uint64_t)p[1]; t2 = (uint64_t)p[2]; - t3 = (uint64_t)p[3]; t4 = (uint64_t)p[4]; - } - - r->d[0] = t0; r->d[1] = t1; r->d[2] = t2; r->d[3] = t3; r->d[4] = t4; + uint64_t w[8]; + mul_wide(w, a->d, a->d); + reduce_wide(r, w); } -#else /* Portable fallback without __int128 */ +#else /* Portable fallback */ -/* Split 64x64 multiply into 32-bit pieces */ static inline void mul64(uint64_t *hi, uint64_t *lo, uint64_t a, uint64_t b) { - uint64_t a_lo = a & 0xFFFFFFFF; - uint64_t a_hi = a >> 32; - uint64_t b_lo = b & 0xFFFFFFFF; - uint64_t b_hi = b >> 32; - - uint64_t ll = a_lo * b_lo; - uint64_t lh = a_lo * b_hi; - uint64_t hl = a_hi * b_lo; - uint64_t hh = a_hi * b_hi; - + uint64_t a_lo = a & 0xFFFFFFFF, a_hi = a >> 32; + uint64_t b_lo = b & 0xFFFFFFFF, b_hi = b >> 32; + uint64_t ll = a_lo * b_lo, lh = a_lo * b_hi, hl = a_hi * b_lo, hh = a_hi * b_hi; uint64_t mid = (ll >> 32) + (lh & 0xFFFFFFFF) + (hl & 0xFFFFFFFF); *lo = (ll & 0xFFFFFFFF) | (mid << 32); *hi = hh + (lh >> 32) + (hl >> 32) + (mid >> 32); } void fe_mul(secp256k1_fe *r, const secp256k1_fe *a, const secp256k1_fe *b) { - /* Portable schoolbook with manual carry tracking */ - const uint64_t M = FE_LIMB_MASK; - /* 2^260 mod p: each limb is 52 bits, so position 5 is at 260 bits. - * 2^260 mod p = 2^4 * (2^256 mod p) = 16 * 0x1000003D1 = 0x10000003D10 */ - const uint64_t R = 0x1000003D10ULL; /* 2^260 mod p = 16 * (2^32 + 977) */ - uint64_t c_hi, c_lo, tmp_hi, tmp_lo; - uint64_t t[5] = {0}; - int i, j; - - /* Simplified portable version - accumulate products */ - for (i = 0; i < 5; i++) { - uint64_t acc_lo = 0, acc_hi = 0; - for (j = 0; j <= i; j++) { - mul64(&tmp_hi, &tmp_lo, a->d[j], b->d[i - j]); - acc_lo += tmp_lo; - acc_hi += tmp_hi + (acc_lo < tmp_lo ? 1 : 0); - } - /* Folded products (j+k >= 5 contribute with factor R) */ - for (j = i + 1; j < 5; j++) { - int k = 5 + i - j; - if (k < 5) { - mul64(&tmp_hi, &tmp_lo, a->d[j], b->d[k]); - /* Multiply by R and add */ - mul64(&c_hi, &c_lo, tmp_lo, R); - acc_lo += c_lo; - acc_hi += c_hi + (acc_lo < c_lo ? 1 : 0); - } - } - t[i] = acc_lo & M; - /* Carry to next limb */ - if (i < 4) { - /* Shift right by 52 */ - uint64_t carry = (acc_lo >> 52) | (acc_hi << 12); - t[i + 1] = carry; + uint64_t w[8] = {0}; + for (int i = 0; i < 4; i++) { + uint64_t carry = 0; + for (int j = 0; j < 4; j++) { + uint64_t hi, lo; + mul64(&hi, &lo, a->d[i], b->d[j]); + uint64_t sum = w[i+j] + lo + carry; + carry = hi + (sum < w[i+j] ? 1 : 0); + w[i+j] = sum; } + w[i+4] += carry; + } + /* Reduce: w[0..3] + w[4..7] * C */ + uint64_t c_lo, c_hi; + uint64_t carry2 = 0; + for (int i = 0; i < 4; i++) { + mul64(&c_hi, &c_lo, w[i+4], FIELD_C); + uint64_t sum = w[i] + c_lo + carry2; + carry2 = c_hi + (sum < w[i] ? 1 : 0); + r->d[i] = sum; + } + if (carry2) { + mul64(&c_hi, &c_lo, carry2, FIELD_C); + uint64_t sum = r->d[0] + c_lo; + r->d[0] = sum; + if (sum < c_lo) { r->d[1]++; if (!r->d[1]) { r->d[2]++; if (!r->d[2]) r->d[3]++; } } } - - r->d[0] = t[0]; r->d[1] = t[1]; r->d[2] = t[2]; r->d[3] = t[3]; r->d[4] = t[4]; fe_normalize(r); } -void fe_sqr(secp256k1_fe *r, const secp256k1_fe *a) { - fe_mul(r, a, a); -} +void fe_sqr(secp256k1_fe *r, const secp256k1_fe *a) { fe_mul(r, a, a); } #endif /* HAVE_INT128 */ @@ -222,162 +184,73 @@ void fe_sqr(secp256k1_fe *r, const secp256k1_fe *a) { static void fe_sqr_n(secp256k1_fe *r, const secp256k1_fe *a, int n) { *r = *a; - for (int i = 0; i < n; i++) { - fe_sqr(r, r); - } + for (int i = 0; i < n; i++) fe_sqr(r, r); } /* ==================== Inversion (Fermat: a^(p-2)) ==================== */ void fe_inv(secp256k1_fe *r, const secp256k1_fe *a) { secp256k1_fe x2, x3, x6, x9, x11, x22, x44, x88, x176, x220, x223; - - fe_sqr(&x2, a); - fe_mul(&x2, &x2, a); - - fe_sqr(&x3, &x2); - fe_mul(&x3, &x3, a); - - fe_sqr_n(&x6, &x3, 3); - fe_mul(&x6, &x6, &x3); - - fe_sqr_n(&x9, &x6, 3); - fe_mul(&x9, &x9, &x3); - - fe_sqr_n(&x11, &x9, 2); - fe_mul(&x11, &x11, &x2); - - fe_sqr_n(&x22, &x11, 11); - fe_mul(&x22, &x22, &x11); - - fe_sqr_n(&x44, &x22, 22); - fe_mul(&x44, &x44, &x22); - - fe_sqr_n(&x88, &x44, 44); - fe_mul(&x88, &x88, &x44); - - fe_sqr_n(&x176, &x88, 88); - fe_mul(&x176, &x176, &x88); - - fe_sqr_n(&x220, &x176, 44); - fe_mul(&x220, &x220, &x44); - - fe_sqr_n(&x223, &x220, 3); - fe_mul(&x223, &x223, &x3); - - fe_sqr_n(r, &x223, 23); - fe_mul(r, r, &x22); - fe_sqr_n(r, r, 5); - fe_mul(r, r, a); - fe_sqr_n(r, r, 3); - fe_mul(r, r, &x2); - fe_sqr_n(r, r, 2); - fe_mul(r, r, a); + fe_sqr(&x2, a); fe_mul(&x2, &x2, a); + fe_sqr(&x3, &x2); fe_mul(&x3, &x3, a); + fe_sqr_n(&x6, &x3, 3); fe_mul(&x6, &x6, &x3); + fe_sqr_n(&x9, &x6, 3); fe_mul(&x9, &x9, &x3); + fe_sqr_n(&x11, &x9, 2); fe_mul(&x11, &x11, &x2); + fe_sqr_n(&x22, &x11, 11); fe_mul(&x22, &x22, &x11); + fe_sqr_n(&x44, &x22, 22); fe_mul(&x44, &x44, &x22); + fe_sqr_n(&x88, &x44, 44); fe_mul(&x88, &x88, &x44); + fe_sqr_n(&x176, &x88, 88); fe_mul(&x176, &x176, &x88); + fe_sqr_n(&x220, &x176, 44); fe_mul(&x220, &x220, &x44); + fe_sqr_n(&x223, &x220, 3); fe_mul(&x223, &x223, &x3); + fe_sqr_n(r, &x223, 23); fe_mul(r, r, &x22); + fe_sqr_n(r, r, 5); fe_mul(r, r, a); + fe_sqr_n(r, r, 3); fe_mul(r, r, &x2); + fe_sqr_n(r, r, 2); fe_mul(r, r, a); } /* ==================== Square root ==================== */ int fe_sqrt(secp256k1_fe *r, const secp256k1_fe *a) { - secp256k1_fe x2, x3, x6, x9, x11, x22, x44, x88, x176, x220, x223; - secp256k1_fe t, check; - - fe_sqr(&x2, a); - fe_mul(&x2, &x2, a); - - fe_sqr(&x3, &x2); - fe_mul(&x3, &x3, a); - - fe_sqr_n(&x6, &x3, 3); - fe_mul(&x6, &x6, &x3); - - fe_sqr_n(&x9, &x6, 3); - fe_mul(&x9, &x9, &x3); - - fe_sqr_n(&x11, &x9, 2); - fe_mul(&x11, &x11, &x2); - - fe_sqr_n(&x22, &x11, 11); - fe_mul(&x22, &x22, &x11); - - fe_sqr_n(&x44, &x22, 22); - fe_mul(&x44, &x44, &x22); - - fe_sqr_n(&x88, &x44, 44); - fe_mul(&x88, &x88, &x44); - - fe_sqr_n(&x176, &x88, 88); - fe_mul(&x176, &x176, &x88); - - fe_sqr_n(&x220, &x176, 44); - fe_mul(&x220, &x220, &x44); - - fe_sqr_n(&x223, &x220, 3); - fe_mul(&x223, &x223, &x3); - - /* (p+1)/4 exponent: same chain but different tail */ - fe_sqr_n(r, &x223, 23); - fe_mul(r, r, &x22); - fe_sqr_n(r, r, 6); - fe_mul(r, r, &x2); + secp256k1_fe x2, x3, x6, x9, x11, x22, x44, x88, x176, x220, x223, check; + fe_sqr(&x2, a); fe_mul(&x2, &x2, a); + fe_sqr(&x3, &x2); fe_mul(&x3, &x3, a); + fe_sqr_n(&x6, &x3, 3); fe_mul(&x6, &x6, &x3); + fe_sqr_n(&x9, &x6, 3); fe_mul(&x9, &x9, &x3); + fe_sqr_n(&x11, &x9, 2); fe_mul(&x11, &x11, &x2); + fe_sqr_n(&x22, &x11, 11); fe_mul(&x22, &x22, &x11); + fe_sqr_n(&x44, &x22, 22); fe_mul(&x44, &x44, &x22); + fe_sqr_n(&x88, &x44, 44); fe_mul(&x88, &x88, &x44); + fe_sqr_n(&x176, &x88, 88); fe_mul(&x176, &x176, &x88); + fe_sqr_n(&x220, &x176, 44); fe_mul(&x220, &x220, &x44); + fe_sqr_n(&x223, &x220, 3); fe_mul(&x223, &x223, &x3); + fe_sqr_n(r, &x223, 23); fe_mul(r, r, &x22); + fe_sqr_n(r, r, 6); fe_mul(r, r, &x2); fe_sqr_n(r, r, 2); - - /* Verify: r^2 == a */ fe_sqr(&check, r); - fe_normalize_full(&check); - t = *a; - fe_normalize_full(&t); - return fe_equal(&check, &t); + return fe_equal(&check, a); } /* ==================== Half ==================== */ void fe_half(secp256k1_fe *r, const secp256k1_fe *a) { - /* - * Compute a/2 mod p. - * If a is even, just shift right by 1. - * If a is odd, add p (which is odd, so a+p is even), then shift right by 1. - * - * We work on the full 256-bit value to avoid carry issues with 5x52 limbs. - * Convert to 4x64, do the conditional add + shift, convert back. - */ secp256k1_fe t = *a; fe_normalize_full(&t); - - /* Reconstruct 4x64 from 5x52 */ - uint64_t v[4]; - v[0] = t.d[0] | (t.d[1] << 52); - v[1] = (t.d[1] >> 12) | (t.d[2] << 40); - v[2] = (t.d[2] >> 24) | (t.d[3] << 28); - v[3] = (t.d[3] >> 36) | (t.d[4] << 16); - - /* p in 4x64 little-endian */ static const uint64_t P[4] = { 0xFFFFFFFEFFFFFC2FULL, 0xFFFFFFFFFFFFFFFFULL, 0xFFFFFFFFFFFFFFFFULL, 0xFFFFFFFFFFFFFFFFULL }; - uint64_t carry = 0; - if (v[0] & 1) { - /* Add p */ + if (t.d[0] & 1) { for (int i = 0; i < 4; i++) { - uint64_t sum = v[i] + P[i] + carry; - carry = (sum < v[i]) || (carry && sum == v[i]) ? 1 : 0; - v[i] = sum; + uint64_t sum = t.d[i] + P[i] + carry; + carry = (sum < t.d[i]) || (carry && sum == t.d[i]) ? 1 : 0; + t.d[i] = sum; } } - - /* Shift right by 1, including the carry bit */ - v[0] = (v[0] >> 1) | (v[1] << 63); - v[1] = (v[1] >> 1) | (v[2] << 63); - v[2] = (v[2] >> 1) | (v[3] << 63); - v[3] = (v[3] >> 1) | (carry << 63); - - /* Convert back to 5x52 */ - r->d[0] = v[0] & FE_LIMB_MASK; - r->d[1] = ((v[0] >> 52) | (v[1] << 12)) & FE_LIMB_MASK; - r->d[2] = ((v[1] >> 40) | (v[2] << 24)) & FE_LIMB_MASK; - r->d[3] = ((v[2] >> 28) | (v[3] << 36)) & FE_LIMB_MASK; - r->d[4] = v[3] >> 16; + r->d[0] = (t.d[0] >> 1) | (t.d[1] << 63); + r->d[1] = (t.d[1] >> 1) | (t.d[2] << 63); + r->d[2] = (t.d[2] >> 1) | (t.d[3] << 63); + r->d[3] = (t.d[3] >> 1) | (carry << 63); } /* ==================== Serialization ==================== */ @@ -385,53 +258,27 @@ void fe_half(secp256k1_fe *r, const secp256k1_fe *a) { void fe_to_bytes(uint8_t *out32, const secp256k1_fe *a) { secp256k1_fe t = *a; fe_normalize_full(&t); - - /* Reconstruct the 256-bit value from 5x52 limbs (little-endian) */ - /* and serialize as big-endian bytes */ - uint64_t v[4]; - v[0] = t.d[0] | (t.d[1] << 52); /* bits 0..103 */ - v[1] = (t.d[1] >> 12) | (t.d[2] << 40); /* bits 64..167 */ - v[2] = (t.d[2] >> 24) | (t.d[3] << 28); /* bits 128..231 */ - v[3] = (t.d[3] >> 36) | (t.d[4] << 16); /* bits 192..255 */ - - /* Write as big-endian */ - for (int i = 0; i < 8; i++) { - out32[31 - i] = (uint8_t)(v[0] >> (i * 8)); - out32[23 - i] = (uint8_t)(v[1] >> (i * 8)); - out32[15 - i] = (uint8_t)(v[2] >> (i * 8)); - out32[7 - i] = (uint8_t)(v[3] >> (i * 8)); + for (int i = 0; i < 4; i++) { + uint64_t v = t.d[3 - i]; + for (int j = 0; j < 8; j++) + out32[i * 8 + j] = (uint8_t)(v >> ((7 - j) * 8)); } } int fe_from_bytes(secp256k1_fe *r, const uint8_t *in32) { - /* Read 32 bytes big-endian into 4x64-bit, then split into 5x52 */ - uint64_t v[4] = {0}; - for (int i = 0; i < 8; i++) { - v[3] |= (uint64_t)in32[i] << ((7 - i) * 8); - v[2] |= (uint64_t)in32[8 + i] << ((7 - i) * 8); - v[1] |= (uint64_t)in32[16 + i] << ((7 - i) * 8); - v[0] |= (uint64_t)in32[24 + i] << ((7 - i) * 8); + for (int i = 0; i < 4; i++) { + uint64_t v = 0; + for (int j = 0; j < 8; j++) + v = (v << 8) | in32[i * 8 + j]; + r->d[3 - i] = v; } - - /* Split 4x64 into 5x52 */ - r->d[0] = v[0] & FE_LIMB_MASK; - r->d[1] = ((v[0] >> 52) | (v[1] << 12)) & FE_LIMB_MASK; - r->d[2] = ((v[1] >> 40) | (v[2] << 24)) & FE_LIMB_MASK; - r->d[3] = ((v[2] >> 28) | (v[3] << 36)) & FE_LIMB_MASK; - r->d[4] = v[3] >> 16; - - /* Check < p */ - secp256k1_fe t = *r; - fe_normalize_full(&t); - /* If normalization changed it, original was >= p */ return 1; } int fe_cmp(const secp256k1_fe *a, const secp256k1_fe *b) { secp256k1_fe ta = *a, tb = *b; - fe_normalize_full(&ta); - fe_normalize_full(&tb); - for (int i = 4; i >= 0; i--) { + fe_normalize_full(&ta); fe_normalize_full(&tb); + for (int i = 3; i >= 0; i--) { if (ta.d[i] < tb.d[i]) return -1; if (ta.d[i] > tb.d[i]) return 1; } diff --git a/quartz/src/main/c/secp256k1/field.h b/quartz/src/main/c/secp256k1/field.h index 0d9c7a0fb2..e429897b08 100644 --- a/quartz/src/main/c/secp256k1/field.h +++ b/quartz/src/main/c/secp256k1/field.h @@ -1,160 +1,135 @@ /* * Copyright (c) 2025 Vitor Pamplona * - * Field arithmetic modulo p = 2^256 - 2^32 - 977 using 5x52-bit limbs. + * Field arithmetic modulo p = 2^256 - 2^32 - 977 using 4x64-bit limbs. * - * Each limb holds up to 52 bits with 12 bits of headroom, allowing - * multiple additions without reduction (lazy reduction). This is the - * key advantage over the Kotlin 4x64-bit representation which requires - * reduction after every add/sub. - * - * On ARM64: uses UMULH/MUL instructions via __int128 - * On x86_64: uses MULQ via __int128 - * Fallback: portable 64-bit C + * Same representation as the Kotlin implementation (Fe4): 4 fully-packed + * 64-bit limbs in little-endian order. This was chosen over the 5x52-bit + * representation because benchmark testing showed fewer multiplies (16 vs 25) + * outweighs the lazy reduction advantage of 5x52 on both JVM and native. */ #ifndef SECP256K1_FIELD_H #define SECP256K1_FIELD_H #include "secp256k1_c.h" -#define FE_LIMB_BITS 52 -#define FE_LIMB_MASK ((uint64_t)0xFFFFFFFFFFFFF) /* 52-bit mask */ +/* ==================== Field Element (4x64-bit limbs, little-endian) ==================== */ -/* ==================== Constants ==================== */ - -static const secp256k1_fe FE_ZERO = {{0, 0, 0, 0, 0}}; -static const secp256k1_fe FE_ONE = {{1, 0, 0, 0, 0}}; - -/* p = 2^256 - 2^32 - 977 in 5x52 limbs */ +/* p = 2^256 - 2^32 - 977 in 4x64 little-endian */ static const secp256k1_fe FE_P = {{ - 0xFFFFEFFFFFC2FULL, /* 4503595332402223 */ - 0xFFFFFFFFFFFFFULL, /* 4503599627370495 */ - 0xFFFFFFFFFFFFFULL, /* 4503599627370495 */ - 0xFFFFFFFFFFFFFULL, /* 4503599627370495 */ - 0x0FFFFFFFFFFFFULL /* 281474976710655 (48-bit top limb) */ + 0xFFFFFFFEFFFFFC2FULL, + 0xFFFFFFFFFFFFFFFFULL, + 0xFFFFFFFFFFFFFFFFULL, + 0xFFFFFFFFFFFFFFFFULL }}; -/* ==================== Core Operations ==================== */ +static const secp256k1_fe FE_ZERO = {{0, 0, 0, 0}}; +static const secp256k1_fe FE_ONE = {{1, 0, 0, 0}}; -/* Normalize to canonical form [0, p) */ -static inline void fe_normalize(secp256k1_fe *r) { - uint64_t t0 = r->d[0], t1 = r->d[1], t2 = r->d[2], t3 = r->d[3], t4 = r->d[4]; - uint64_t m; +/* P[0] cached for hot path */ +#define FE_P0 0xFFFFFFFEFFFFFC2FULL - /* Reduce carries */ - t1 += t0 >> 52; t0 &= FE_LIMB_MASK; - t2 += t1 >> 52; t1 &= FE_LIMB_MASK; - t3 += t2 >> 52; t2 &= FE_LIMB_MASK; - t4 += t3 >> 52; t3 &= FE_LIMB_MASK; - - /* t4 may overflow 48 bits; fold top bits: 2^256 = 2^32 + 977 (mod p) */ - m = t4 >> 48; - t4 &= 0xFFFFFFFFFFFFULL; /* 48-bit mask */ - t0 += m * 0x1000003D1ULL; - t1 += t0 >> 52; t0 &= FE_LIMB_MASK; - t2 += t1 >> 52; t1 &= FE_LIMB_MASK; - t3 += t2 >> 52; t2 &= FE_LIMB_MASK; - t4 += t3 >> 52; t3 &= FE_LIMB_MASK; - - /* Final conditional subtraction of p */ - /* p in 5x52: [0xFFFFEFFFFFC2F, 0xFFFFFFFFFFFFF, 0xFFFFFFFFFFFFF, 0xFFFFFFFFFFFFF, 0x0FFFFFFFFFFFF] */ - m = (t4 == 0x0FFFFFFFFFFFFULL) & - (t3 == FE_LIMB_MASK) & - (t2 == FE_LIMB_MASK) & - (t1 == FE_LIMB_MASK) & - (t0 >= 0xFFFFEFFFFFC2FULL); - t0 -= m * 0xFFFFEFFFFFC2FULL; - t1 -= m * FE_LIMB_MASK; - t2 -= m * FE_LIMB_MASK; - t3 -= m * FE_LIMB_MASK; - t4 -= m * 0x0FFFFFFFFFFFFULL; - - /* Re-propagate borrows */ - if (m) { - /* After subtracting p, no borrows are possible if t >= p */ - /* But handle just in case of numerical edge cases */ - } - - r->d[0] = t0; r->d[1] = t1; r->d[2] = t2; r->d[3] = t3; r->d[4] = t4; -} - -/* Normalize fully (for comparison/serialization) */ -static inline void fe_normalize_full(secp256k1_fe *r) { - fe_normalize(r); - fe_normalize(r); /* Second pass for edge cases */ -} +/* ==================== Inline Helpers ==================== */ static inline int fe_is_zero(const secp256k1_fe *a) { + /* Normalize before checking — elements may be unreduced */ secp256k1_fe t = *a; - fe_normalize_full(&t); - return (t.d[0] | t.d[1] | t.d[2] | t.d[3] | t.d[4]) == 0; + /* Quick check: if all limbs are in range and < p, it's normalized */ + if (t.d[3] == UINT64_MAX && t.d[2] == UINT64_MAX && + t.d[1] == UINT64_MAX && t.d[0] >= FE_P0) { + /* >= p, reduce */ + t.d[0] -= FE_P0; t.d[1] = 0; t.d[2] = 0; t.d[3] = 0; + } + return (t.d[0] | t.d[1] | t.d[2] | t.d[3]) == 0; } static inline int fe_equal(const secp256k1_fe *a, const secp256k1_fe *b) { secp256k1_fe ta = *a, tb = *b; - fe_normalize_full(&ta); - fe_normalize_full(&tb); - return (ta.d[0] == tb.d[0]) & (ta.d[1] == tb.d[1]) & (ta.d[2] == tb.d[2]) & - (ta.d[3] == tb.d[3]) & (ta.d[4] == tb.d[4]); + /* Normalize both */ + if (ta.d[3] == UINT64_MAX && ta.d[2] == UINT64_MAX && + ta.d[1] == UINT64_MAX && ta.d[0] >= FE_P0) { + ta.d[0] -= FE_P0; ta.d[1] = 0; ta.d[2] = 0; ta.d[3] = 0; + } + if (tb.d[3] == UINT64_MAX && tb.d[2] == UINT64_MAX && + tb.d[1] == UINT64_MAX && tb.d[0] >= FE_P0) { + tb.d[0] -= FE_P0; tb.d[1] = 0; tb.d[2] = 0; tb.d[3] = 0; + } + return (ta.d[0] == tb.d[0]) & (ta.d[1] == tb.d[1]) & + (ta.d[2] == tb.d[2]) & (ta.d[3] == tb.d[3]); } static inline int fe_is_odd(const secp256k1_fe *a) { secp256k1_fe t = *a; - fe_normalize_full(&t); + if (t.d[3] == UINT64_MAX && t.d[2] == UINT64_MAX && + t.d[1] == UINT64_MAX && t.d[0] >= FE_P0) { + t.d[0] -= FE_P0; t.d[1] = 0; t.d[2] = 0; t.d[3] = 0; + } return (int)(t.d[0] & 1); } -/* r = a + b (lazy: no reduction) */ +/* Normalize: if a >= p, subtract p. Inline for hot path. */ +static inline void fe_normalize(secp256k1_fe *a) { + if (a->d[3] == UINT64_MAX && a->d[2] == UINT64_MAX && + a->d[1] == UINT64_MAX && a->d[0] >= FE_P0) { + a->d[0] -= FE_P0; + a->d[1] = 0; + a->d[2] = 0; + a->d[3] = 0; + } +} + +static inline void fe_normalize_full(secp256k1_fe *a) { + fe_normalize(a); +} + +/* r = a + b mod p */ static inline void fe_add(secp256k1_fe *r, const secp256k1_fe *a, const secp256k1_fe *b) { - r->d[0] = a->d[0] + b->d[0]; - r->d[1] = a->d[1] + b->d[1]; - r->d[2] = a->d[2] + b->d[2]; - r->d[3] = a->d[3] + b->d[3]; - r->d[4] = a->d[4] + b->d[4]; + uint64_t carry = 0; + for (int i = 0; i < 4; i++) { + uint64_t sum = a->d[i] + b->d[i] + carry; + carry = (sum < a->d[i]) || (carry && sum == a->d[i]) ? 1 : 0; + r->d[i] = sum; + } + if (carry) { + /* Overflow past 2^256: add 2^256 mod p = C = 0x1000003D1 */ + uint64_t s = r->d[0] + 0x1000003D1ULL; + uint64_t c = (s < r->d[0]) ? 1 : 0; + r->d[0] = s; + if (c) { r->d[1]++; if (!r->d[1]) { r->d[2]++; if (!r->d[2]) r->d[3]++; } } + } + fe_normalize(r); } -/* r += a (in-place lazy add) */ +/* r += a */ static inline void fe_add_assign(secp256k1_fe *r, const secp256k1_fe *a) { - r->d[0] += a->d[0]; - r->d[1] += a->d[1]; - r->d[2] += a->d[2]; - r->d[3] += a->d[3]; - r->d[4] += a->d[4]; + secp256k1_fe t = *r; + fe_add(r, &t, a); } -/* r = -a mod p. Computes (m+1)*p - a to keep limbs positive (works for magnitude <= m) */ +/* r = -a mod p = P - a */ static inline void fe_negate(secp256k1_fe *r, const secp256k1_fe *a, int m) { - /* Add (m+1)*p and subtract a */ - uint64_t mp = (uint64_t)(m + 1); - r->d[0] = mp * 0xFFFFEFFFFFC2FULL - a->d[0]; - r->d[1] = mp * 0xFFFFFFFFFFFFFULL - a->d[1]; - r->d[2] = mp * 0xFFFFFFFFFFFFFULL - a->d[2]; - r->d[3] = mp * 0xFFFFFFFFFFFFFULL - a->d[3]; - r->d[4] = mp * 0x0FFFFFFFFFFFFULL - a->d[4]; + (void)m; /* magnitude parameter not needed for 4x64 */ + if (fe_is_zero(a)) { + *r = FE_ZERO; + return; + } + uint64_t borrow = 0; + for (int i = 0; i < 4; i++) { + uint64_t diff = FE_P.d[i] - a->d[i] - borrow; + borrow = (FE_P.d[i] < a->d[i] + borrow) || (borrow && a->d[i] == UINT64_MAX) ? 1 : 0; + r->d[i] = diff; + } } -/* r = a * b mod p */ +/* ==================== Function declarations ==================== */ + void fe_mul(secp256k1_fe *r, const secp256k1_fe *a, const secp256k1_fe *b); - -/* r = a^2 mod p */ void fe_sqr(secp256k1_fe *r, const secp256k1_fe *a); - -/* r = a^(-1) mod p (Fermat: a^(p-2)) */ void fe_inv(secp256k1_fe *r, const secp256k1_fe *a); - -/* r = sqrt(a) mod p, returns 1 on success */ int fe_sqrt(secp256k1_fe *r, const secp256k1_fe *a); - -/* r = a/2 mod p */ void fe_half(secp256k1_fe *r, const secp256k1_fe *a); - -/* Serialize field element to 32-byte big-endian */ void fe_to_bytes(uint8_t *out32, const secp256k1_fe *a); - -/* Deserialize 32-byte big-endian to field element */ int fe_from_bytes(secp256k1_fe *r, const uint8_t *in32); - -/* Compare field elements: -1, 0, 1 */ int fe_cmp(const secp256k1_fe *a, const secp256k1_fe *b); #endif /* SECP256K1_FIELD_H */ diff --git a/quartz/src/main/c/secp256k1/point.c b/quartz/src/main/c/secp256k1/point.c index 0ca4851eab..7379e7f8d8 100644 --- a/quartz/src/main/c/secp256k1/point.c +++ b/quartz/src/main/c/secp256k1/point.c @@ -11,21 +11,23 @@ /* G_x = 0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798 */ /* G_y = 0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8 */ +/* 4x64 little-endian */ const secp256k1_ge SECP256K1_G = { - .x = {{0x2815B16F81798ULL, 0xDB2DCE28D959FULL, 0xE870B07029BFCULL, - 0xBBAC55A06295CULL, 0x079BE667EF9DCULL}}, - .y = {{0x7D08FFB10D4B8ULL, 0x48A68554199C4ULL, 0xE1108A8FD17B4ULL, - 0xC4655DA4FBFC0ULL, 0x0483ADA7726A3ULL}} + .x = {{0x59F2815B16F81798ULL, 0x029BFCDB2DCE28D9ULL, + 0x55A06295CE870B07ULL, 0x79BE667EF9DCBBACULL}}, + .y = {{0x9C47D08FFB10D4B8ULL, 0xFD17B448A6855419ULL, + 0x5DA4FBFC0E1108A8ULL, 0x483ADA7726A3C465ULL}} }; -/* GLV beta: cube root of unity mod p (5x52 limbs) */ +/* GLV beta: cube root of unity mod p (4x64 little-endian) */ +/* beta = 0x7AE96A2B657C07106E64479EAC3434E99CF0497512F58995C1396C28719501EE */ static const secp256k1_fe GLV_BETA = {{ - 0x96C28719501EEULL, 0x7512F58995C13ULL, 0xC3434E99CF049ULL, - 0x07106E64479EAULL, 0x07AE96A2B657CULL + 0xC1396C28719501EEULL, 0x9CF0497512F58995ULL, + 0x6E64479EAC3434E9ULL, 0x7AE96A2B657C0710ULL }}; /* Curve constant b = 7 */ -static const secp256k1_fe FE_SEVEN = {{7, 0, 0, 0, 0}}; +static const secp256k1_fe FE_SEVEN = {{7, 0, 0, 0}}; /* ==================== Precomputed Tables ==================== */ @@ -67,6 +69,13 @@ int gej_is_infinity(const secp256k1_gej *r) { void gej_double(secp256k1_gej *r, const secp256k1_gej *p) { secp256k1_fe s, l, t, u; + /* Handle aliasing: if r == p, copy input first */ + secp256k1_gej tmp; + if (r == p) { + tmp = *p; + p = &tmp; + } + if (p->infinity) { gej_set_infinity(r); return; diff --git a/quartz/src/main/c/secp256k1/secp256k1_c.h b/quartz/src/main/c/secp256k1/secp256k1_c.h index 7d42e80377..b60b46825a 100644 --- a/quartz/src/main/c/secp256k1/secp256k1_c.h +++ b/quartz/src/main/c/secp256k1/secp256k1_c.h @@ -65,21 +65,13 @@ /* * Field element modulo p = 2^256 - 2^32 - 977. - * 5 limbs of 52 bits each, with 12 bits of headroom per limb. - * This allows 3-8 chained additions without reduction (lazy reduction). - * - * Magnitude tracking: after N additions without reduction, each limb - * can be up to N * 2^52. We reduce when magnitude exceeds safe limits. + * 4 limbs of 64 bits each, fully packed, little-endian. + * Same representation as the Kotlin Fe4 class. */ typedef struct { - uint64_t d[5]; + uint64_t d[4]; } secp256k1_fe; -/* Wide result of field multiplication (used internally) */ -typedef struct { - uint64_t d[10]; -} secp256k1_fe_wide; - /* ==================== Scalar (mod n) ==================== */ typedef struct {