From 1dd8edcd414a238efbd529b536acc6bf9e9ccea0 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 7 Apr 2026 23:14:11 +0000 Subject: [PATCH] perf: optimize secp256k1 field arithmetic and point operations for Android MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Unroll hot-path loops in U256 (mulWide, sqrWide, addTo, subTo) to eliminate loop control overhead and array bounds checks that ART JIT on Android does not optimize as aggressively as HotSpot. These functions are called ~1,900× per signature verification. Optimize unsignedMultiplyHighFallback to compute the unsigned 128-bit product directly from 32-bit sub-products, avoiding the signed multiplyHigh + correction path. Saves ~8 instructions per call on Android < API 31 (~30,000 calls per verify). Unroll FieldP.reduceWide carry propagation and inline the reduceSelf subtraction (avoiding P array access and U256.subTo call). Unroll FieldP.half with P[1..3]=-1 inlined as mask. Replace copyFrom pattern in mul, mulG, and mulDoubleG with ping-pong point buffers — alternate between two MutablePoints via reference swap instead of copying 12 Longs after every addition. Eliminates ~170 copyFrom calls per verify and avoids the internal copy buffer in doublePoint's out===inp path (~130 per verify). Also removes the MutablePoint() allocation in mulG (3 LongArrays per sign/pubCreate). https://claude.ai/code/session_017UbWduFi1sLUsgVUMUH2nx --- .../quartz/utils/secp256k1/FieldP.kt | 190 +++++--- .../quartz/utils/secp256k1/MultiplyHigh.kt | 20 +- .../quartz/utils/secp256k1/Point.kt | 112 +++-- .../quartz/utils/secp256k1/U256.kt | 438 ++++++++++++++---- 4 files changed, 586 insertions(+), 174 deletions(-) diff --git a/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/FieldP.kt b/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/FieldP.kt index 7fdfdc0af..555b3b3e8 100644 --- a/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/FieldP.kt +++ b/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/FieldP.kt @@ -67,14 +67,14 @@ internal object FieldP { val carry = U256.addTo(out, a, b) if (carry != 0) { // Overflow past 2^256: add 2^256 mod p = 2^32 + 977 = 0x1000003D1 - // This fits in 33 bits. Add to limb[0] with carry propagation. - val s1 = out[0] + 4294968273L // 2^32 + 977 + val s1 = out[0] + 4294968273L val c1 = if (s1.toULong() < out[0].toULong()) 1L else 0L out[0] = s1 if (c1 != 0L) { - for (i in 1 until 4) { - out[i]++ - if (out[i] != 0L) break + out[1]++ + if (out[1] == 0L) { + out[2]++ + if (out[2] == 0L) out[3]++ } } } @@ -145,26 +145,51 @@ internal object FieldP { /** * out = a / 2 mod p. Branchless: if odd, add p first (p is odd → a+p is even). + * Unrolled, with P[1..3]=-1 inlined as `mask` (since -1 & mask = mask). */ fun half( out: LongArray, a: LongArray, ) { val mask = -(a[0] and 1L) // all 1s if odd, all 0s if even - var carry = 0L - for (i in 0 until 4) { - val pMasked = P[i] and mask - val s1 = a[i] + pMasked - val c1 = if (s1.toULong() < a[i].toULong()) 1L else 0L - val s2 = s1 + carry - val c2 = if (s2.toULong() < s1.toULong()) 1L else 0L - out[i] = s2 - carry = c1 + c2 - } - // Right-shift by 1 - for (i in 0 until 3) { - out[i] = (out[i] ushr 1) or (out[i + 1] shl 63) - } + val p0 = P0 and mask // P[0] masked; P[1..3] are -1, so P[i]&mask = mask + var s1: Long + var s2: Long + var c1: Long + var c2: Long + + // Conditional add: out = a + (P & mask), unrolled + // Limb 0 + s1 = a[0] + p0 + c1 = if (s1.toULong() < a[0].toULong()) 1L else 0L + out[0] = s1 + var carry = c1 + // Limb 1 + s1 = a[1] + mask + c1 = if (s1.toULong() < a[1].toULong()) 1L else 0L + s2 = s1 + carry + c2 = if (s2.toULong() < s1.toULong()) 1L else 0L + out[1] = s2 + carry = c1 + c2 + // Limb 2 + s1 = a[2] + mask + c1 = if (s1.toULong() < a[2].toULong()) 1L else 0L + s2 = s1 + carry + c2 = if (s2.toULong() < s1.toULong()) 1L else 0L + out[2] = s2 + carry = c1 + c2 + // Limb 3 + s1 = a[3] + mask + c1 = if (s1.toULong() < a[3].toULong()) 1L else 0L + s2 = s1 + carry + c2 = if (s2.toULong() < s1.toULong()) 1L else 0L + out[3] = s2 + carry = c1 + c2 + + // Right-shift by 1 (unrolled) + out[0] = (out[0] ushr 1) or (out[1] shl 63) + out[1] = (out[1] ushr 1) or (out[2] shl 63) + out[2] = (out[2] ushr 1) or (out[3] shl 63) out[3] = (out[3] ushr 1) or (carry shl 63) } @@ -298,81 +323,122 @@ internal object FieldP { // ==================== Reduction ==================== + // P[0] cached as a constant to avoid array load in the hot reduceSelf path. + private const val P0 = -4294968273L // 0xFFFFFFFEFFFFFC2F + fun reduceSelf(a: LongArray) { // Exploit P's structure: P = [P0, -1, -1, -1] where P[1..3] = 0xFFFFFFFFFFFFFFFF. // a >= P only if a[3]==a[2]==a[1]==-1 AND a[0] >= P[0]. The first check (a[3]==-1) // fails >99.99% of the time for random field elements, making this a single branch. if (a[3] == -1L && a[2] == -1L && a[1] == -1L && - (a[0] xor Long.MIN_VALUE) >= (P[0] xor Long.MIN_VALUE) + (a[0] xor Long.MIN_VALUE) >= (P0 xor Long.MIN_VALUE) ) { - U256.subTo(a, a, P) + // Inline P subtraction: when a[1..3] = -1 and a[0] >= P0, + // a - P = [a[0] - P0, 0, 0, 0] (no borrows since P[1..3] = -1). + a[0] -= P0 + a[1] = 0L + a[2] = 0L + a[3] = 0L } } /** - * Reduce 512-bit value mod p. + * Reduce 512-bit value mod p. Fully unrolled for ART JIT. * * Uses hi × 2^256 ≡ hi × C (mod p) where C = 2^32 + 977 = 4294968273. - * Since C < 2^33, hi[i] × C fits in 97 bits. We use unsignedMultiplyHigh - * to get the upper 64 bits of each limb×C product. - * - * Three stages: - * 1. Fold 512→~260 bits: lo + hi × C, producing at most ~34-bit carry - * 2. Fold carry × C back into limb[0..3]; propagate carries (may overflow 256 bits) - * 3. If round 2 overflowed, fold the single-bit overflow (≡ C) once more - * Final reduceSelf handles the at-most-one subtraction of p. + * Three stages: fold 512→~260 bits, fold carry×C, final reduceSelf. */ fun reduceWide( out: LongArray, w: LongArray, ) { - // Round 1: acc = lo + hi × C val c = 4294968273L // 2^32 + 977 - var carry = 0L - for (i in 0 until 4) { - val hcLo = w[i + 4] * c - val hcHi = unsignedMultiplyHigh(w[i + 4], c) + var hcLo: Long + var hcHi: Long + var s1: Long + var s2: Long + var c1: Long + var c2: Long - // acc = w[i] + hcLo + carry - val s1 = w[i] + hcLo - val c1 = if (s1.toULong() < w[i].toULong()) 1L else 0L - val s2 = s1 + carry - val c2 = if (s2.toULong() < s1.toULong()) 1L else 0L - out[i] = s2 - carry = hcHi + c1 + c2 - } + // Round 1: acc = lo + hi × C (4 limbs, unrolled) + + // Limb 0 (no carry input) + hcLo = w[4] * c + hcHi = unsignedMultiplyHigh(w[4], c) + s1 = w[0] + hcLo + c1 = if (s1.toULong() < w[0].toULong()) 1L else 0L + out[0] = s1 + var carry = hcHi + c1 + + // Limb 1 + hcLo = w[5] * c + hcHi = unsignedMultiplyHigh(w[5], c) + s1 = w[1] + hcLo + c1 = if (s1.toULong() < w[1].toULong()) 1L else 0L + s2 = s1 + carry + c2 = if (s2.toULong() < s1.toULong()) 1L else 0L + out[1] = s2 + carry = hcHi + c1 + c2 + + // Limb 2 + hcLo = w[6] * c + hcHi = unsignedMultiplyHigh(w[6], c) + s1 = w[2] + hcLo + c1 = if (s1.toULong() < w[2].toULong()) 1L else 0L + s2 = s1 + carry + c2 = if (s2.toULong() < s1.toULong()) 1L else 0L + out[2] = s2 + carry = hcHi + c1 + c2 + + // Limb 3 + hcLo = w[7] * c + hcHi = unsignedMultiplyHigh(w[7], c) + s1 = w[3] + hcLo + c1 = if (s1.toULong() < w[3].toULong()) 1L else 0L + s2 = s1 + carry + c2 = if (s2.toULong() < s1.toULong()) 1L else 0L + out[3] = s2 + carry = hcHi + c1 + c2 // Round 2: if carry > 0, fold carry × C back in if (carry != 0L) { val ccLo = carry * c val ccHi = unsignedMultiplyHigh(carry, c) - val s1 = out[0] + ccLo - val c1 = if (s1.toULong() < out[0].toULong()) 1L else 0L + s1 = out[0] + ccLo + c1 = if (s1.toULong() < out[0].toULong()) 1L else 0L out[0] = s1 + // Propagate carry (unrolled, with early exit) var prop = ccHi + c1 - for (i in 1 until 4) { - if (prop == 0L) break - val s = out[i] + prop - prop = if (s.toULong() < out[i].toULong()) 1L else 0L - out[i] = s - } - // Round 2 carry propagation may overflow past 256 bits. - // This happens when out[0..3] were all 0xFF..FF and the add cascades. - // Overflow of 1 means 2^256 ≡ C (mod p), so add C to out[0..3]. if (prop != 0L) { - val s2 = out[0] + c - val c2 = if (s2.toULong() < out[0].toULong()) 1L else 0L - out[0] = s2 - if (c2 != 0L) { - for (i in 1 until 4) { - out[i]++ - if (out[i] != 0L) break + s1 = out[1] + prop + prop = if (s1.toULong() < out[1].toULong()) 1L else 0L + out[1] = s1 + if (prop != 0L) { + s1 = out[2] + prop + prop = if (s1.toULong() < out[2].toULong()) 1L else 0L + out[2] = s1 + if (prop != 0L) { + s1 = out[3] + prop + prop = if (s1.toULong() < out[3].toULong()) 1L else 0L + out[3] = s1 + } + } + } + // Overflow past 256 bits: 2^256 ≡ C (mod p) + if (prop != 0L) { + s1 = out[0] + c + c1 = if (s1.toULong() < out[0].toULong()) 1L else 0L + out[0] = s1 + if (c1 != 0L) { + out[1]++ + if (out[1] == 0L) { + out[2]++ + if (out[2] == 0L) out[3]++ } } } } - // Final: at most one subtraction of p reduceSelf(out) } diff --git a/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/MultiplyHigh.kt b/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/MultiplyHigh.kt index 4066ec598..77992e782 100644 --- a/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/MultiplyHigh.kt +++ b/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/MultiplyHigh.kt @@ -45,12 +45,28 @@ internal expect fun unsignedMultiplyHigh( ): Long /** - * Fallback: unsigned multiply high from signed multiply high + correction. + * Fallback: unsigned multiply high computed directly from 32-bit sub-products. + * + * Unlike the old approach (signed multiplyHigh + correction), this computes the + * unsigned result directly, avoiding the signed correction branches (if a < 0, + * if b < 0) and the unsigned correction terms (+ (a & (b >> 63)) + (b & (a >> 63))). + * Saves ~8 instructions per call on Android < API 31, where this is the hot path + * (~30,000 calls per signature verify). */ internal fun unsignedMultiplyHighFallback( a: Long, b: Long, -): Long = multiplyHigh(a, b) + (a and (b shr 63)) + (b and (a shr 63)) +): Long { + val aLo = a and 0xFFFFFFFFL + val aHi = a ushr 32 + val bLo = b and 0xFFFFFFFFL + val bHi = b ushr 32 + val mid1 = aHi * bLo + val mid2 = aLo * bHi + val low = aLo * bLo + val carry = ((low ushr 32) + (mid1 and 0xFFFFFFFFL) + (mid2 and 0xFFFFFFFFL)) ushr 32 + return (aHi * bHi) + (mid1 ushr 32) + (mid2 ushr 32) + carry +} /** * Pure-Kotlin fallback for multiplyHigh, using four 32-bit sub-products. diff --git a/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/Point.kt b/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/Point.kt index 427c39006..cb90bd61c 100644 --- a/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/Point.kt +++ b/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/Point.kt @@ -596,15 +596,31 @@ internal object ECPoint { bits-- } - out.setInfinity() - val tmp = s.mixTmp + // Ping-pong: alternate between two point buffers to avoid copyFrom after + // every addition. Saves ~20 copyFroms per call (each = 12 Long copies). + // Also avoids the internal copy in doublePoint (out===inp path). + var cur = out + var alt = s.mixTmp + cur.setInfinity() val negY = s.mixNegY for (i in bits - 1 downTo 0) { - doublePoint(out, out, s) - addWnafMixed(out, tmp, negY, wnaf1, i, pOdd, split.negK1, s) - addWnafMixed(out, tmp, negY, wnaf2, i, pLamOdd, split.negK2, s) + doublePoint(alt, cur, s) + var t = cur + cur = alt + alt = t + if (addWnafMixedPP(cur, alt, negY, wnaf1, i, pOdd, split.negK1, s)) { + t = cur + cur = alt + alt = t + } + if (addWnafMixedPP(cur, alt, negY, wnaf2, i, pLamOdd, split.negK2, s)) { + t = cur + cur = alt + alt = t + } } + if (cur !== out) out.copyFrom(cur) } /** @@ -628,12 +644,20 @@ internal object ECPoint { val s = scratch.get() val table = combTable - out.setInfinity() - val tmp = MutablePoint() + + // Ping-pong: alternate between out and s.mixTmp to avoid copyFrom after + // every addMixed and the internal copy in in-place doublePoint. + // Also eliminates the MutablePoint() allocation that was here before. + var cur = out + var alt = s.mixTmp + cur.setInfinity() for (combOff in COMB_SPACING - 1 downTo 0) { if (combOff < COMB_SPACING - 1) { - doublePoint(out, out, s) + doublePoint(alt, cur, s) + val t = cur + cur = alt + alt = t } for (block in 0 until COMB_BLOCKS) { var mask = 0 @@ -645,11 +669,14 @@ internal object ECPoint { } if (mask != 0) { val entry = table[block * COMB_POINTS + mask] - addMixed(tmp, out, entry.x, entry.y, s) - out.copyFrom(tmp) + addMixed(alt, cur, entry.x, entry.y, s) + val t = cur + cur = alt + alt = t } } } + if (cur !== out) out.copyFrom(cur) } /** @@ -717,48 +744,75 @@ internal object ECPoint { bits-- } - out.setInfinity() - val tmp = sc.mixTmp + // Ping-pong: alternate between out and sc.mixTmp to avoid copyFrom after + // every addition (~170 copies per verify → at most 1). Also avoids the + // internal copy buffer in doublePoint's out===inp path (~130 per verify). + var cur = out + var alt = sc.mixTmp + cur.setInfinity() val negY = sc.mixNegY for (i in bits - 1 downTo 0) { - doublePoint(out, out, sc) + doublePoint(alt, cur, sc) + var t = cur + cur = alt + alt = t // Streams 1-2: G-side (affine tables, mixed addition) - addWnafMixed(out, tmp, negY, wnafS1, i, gOdd, sSplit.negK1, sc) - addWnafMixed(out, tmp, negY, wnafS2, i, gLam, sSplit.negK2, sc) + if (addWnafMixedPP(cur, alt, negY, wnafS1, i, gOdd, sSplit.negK1, sc)) { + t = cur + cur = alt + alt = t + } + if (addWnafMixedPP(cur, alt, negY, wnafS2, i, gLam, sSplit.negK2, sc)) { + t = cur + cur = alt + alt = t + } // Streams 3-4: P-side (affine tables via effective-affine, mixed addition) - addWnafMixed(out, tmp, negY, wnafE1, i, pOdd, eSplit.negK1, sc) - addWnafMixed(out, tmp, negY, wnafE2, i, pLamOdd, eSplit.negK2, sc) + if (addWnafMixedPP(cur, alt, negY, wnafE1, i, pOdd, eSplit.negK1, sc)) { + t = cur + cur = alt + alt = t + } + if (addWnafMixedPP(cur, alt, negY, wnafE2, i, pLamOdd, eSplit.negK2, sc)) { + t = cur + cur = alt + alt = t + } } + if (cur !== out) out.copyFrom(cur) } /** - * Process one wNAF digit with mixed addition. - * The effective sign is: (wNAF digit sign) XOR (GLV negation flag). - * Positive = add as-is, negative = negate the table entry's y. + * Process one wNAF digit with mixed addition (ping-pong version). + * Reads from `cur`, writes result to `alt`. Returns true if an addition was + * performed (caller should swap cur/alt references). + * + * This avoids the copyFrom after every addition — the caller swaps references + * instead (free: just local variable reassignment). */ - private fun addWnafMixed( - out: MutablePoint, - tmp: MutablePoint, + private fun addWnafMixedPP( + cur: MutablePoint, + alt: MutablePoint, negY: LongArray, wnafDigits: IntArray, bitIndex: Int, table: Array, glvNeg: Boolean, s: PointScratch, - ) { - if (bitIndex >= wnafDigits.size) return + ): Boolean { + if (bitIndex >= wnafDigits.size) return false val d = wnafDigits[bitIndex] - if (d == 0) return + if (d == 0) return false val idx = (if (d > 0) d else -d) / 2 val effectiveNeg = (d < 0) xor glvNeg if (!effectiveNeg) { - addMixed(tmp, out, table[idx].x, table[idx].y, s) + addMixed(alt, cur, table[idx].x, table[idx].y, s) } else { FieldP.neg(negY, table[idx].y) - addMixed(tmp, out, table[idx].x, negY, s) + addMixed(alt, cur, table[idx].x, negY, s) } - out.copyFrom(tmp) + return true } /** Process one wNAF digit with full Jacobian addition (for P-side tables). */ diff --git a/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/U256.kt b/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/U256.kt index 421afb462..61a756807 100644 --- a/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/U256.kt +++ b/quartz/src/commonMain/kotlin/com/vitorpamplona/quartz/utils/secp256k1/U256.kt @@ -64,131 +64,407 @@ internal object U256 { return 0 } - /** out = a + b. Returns carry (0 or 1). Safe for aliasing. */ + /** out = a + b. Returns carry (0 or 1). Safe for aliasing. Unrolled for ART JIT. */ fun addTo( out: LongArray, a: LongArray, b: LongArray, ): Int { - var carry = 0L - for (i in 0 until 4) { - val s1 = a[i] + b[i] - val c1 = if (s1.toULong() < a[i].toULong()) 1L else 0L - val s2 = s1 + carry - val c2 = if (s2.toULong() < s1.toULong()) 1L else 0L - out[i] = s2 - carry = c1 + c2 - } + var s1: Long + var s2: Long + var c1: Long + var c2: Long + + // Limb 0 (no carry input) + s1 = a[0] + b[0] + c1 = if (s1.toULong() < a[0].toULong()) 1L else 0L + out[0] = s1 + var carry = c1 + + // Limb 1 + s1 = a[1] + b[1] + c1 = if (s1.toULong() < a[1].toULong()) 1L else 0L + s2 = s1 + carry + c2 = if (s2.toULong() < s1.toULong()) 1L else 0L + out[1] = s2 + carry = c1 + c2 + + // Limb 2 + s1 = a[2] + b[2] + c1 = if (s1.toULong() < a[2].toULong()) 1L else 0L + s2 = s1 + carry + c2 = if (s2.toULong() < s1.toULong()) 1L else 0L + out[2] = s2 + carry = c1 + c2 + + // Limb 3 + s1 = a[3] + b[3] + c1 = if (s1.toULong() < a[3].toULong()) 1L else 0L + s2 = s1 + carry + c2 = if (s2.toULong() < s1.toULong()) 1L else 0L + out[3] = s2 + carry = c1 + c2 + return carry.toInt() } - /** out = a - b. Returns borrow (0 or 1). Safe for aliasing. */ + /** out = a - b. Returns borrow (0 or 1). Safe for aliasing. Unrolled for ART JIT. */ fun subTo( out: LongArray, a: LongArray, b: LongArray, ): Int { - var borrow = 0L - for (i in 0 until 4) { - val d1 = a[i] - b[i] - val c1 = if (a[i].toULong() < b[i].toULong()) 1L else 0L - val d2 = d1 - borrow - val c2 = if (d1.toULong() < borrow.toULong()) 1L else 0L - out[i] = d2 - borrow = c1 + c2 - } + var d1: Long + var d2: Long + var c1: Long + var c2: Long + + // Limb 0 (no borrow input) + d1 = a[0] - b[0] + c1 = if (a[0].toULong() < b[0].toULong()) 1L else 0L + out[0] = d1 + var borrow = c1 + + // Limb 1 + d1 = a[1] - b[1] + c1 = if (a[1].toULong() < b[1].toULong()) 1L else 0L + d2 = d1 - borrow + c2 = if (d1.toULong() < borrow.toULong()) 1L else 0L + out[1] = d2 + borrow = c1 + c2 + + // Limb 2 + d1 = a[2] - b[2] + c1 = if (a[2].toULong() < b[2].toULong()) 1L else 0L + d2 = d1 - borrow + c2 = if (d1.toULong() < borrow.toULong()) 1L else 0L + out[2] = d2 + borrow = c1 + c2 + + // Limb 3 + d1 = a[3] - b[3] + c1 = if (a[3].toULong() < b[3].toULong()) 1L else 0L + d2 = d1 - borrow + c2 = if (d1.toULong() < borrow.toULong()) 1L else 0L + out[3] = d2 + borrow = c1 + c2 + return borrow.toInt() } /** * 4×4 schoolbook multiplication: out = a × b (512-bit result in LongArray(8)). * - * Uses unsignedMultiplyHigh for the upper 64 bits of each 64×64→128-bit product. - * On JVM, this is a hardware intrinsic (single instruction). Total: 16 products - * vs 64 for the previous 8×32-bit representation. + * Fully unrolled: all 16 products are explicit, eliminating loop control overhead + * and array bounds checks. This significantly helps ART JIT on Android, which is + * less aggressive at loop optimization than HotSpot. Called ~1,900× per verify. */ fun mulWide( out: LongArray, a: LongArray, b: LongArray, ) { - for (i in 0 until 8) out[i] = 0L + val a0 = a[0] + val a1 = a[1] + val a2 = a[2] + val a3 = a[3] + val b0 = b[0] + val b1 = b[1] + val b2 = b[2] + val b3 = b[3] + var lo: Long + var hi: Long + var prev: Long + var s: Long + var c1: Long + var c2: Long + var carry: Long - for (i in 0 until 4) { - var carry = 0L - val ai = a[i] - for (j in 0 until 4) { - val lo = ai * b[j] - val hi = unsignedMultiplyHigh(ai, b[j]) + // Row 0: a0 × [b0,b1,b2,b3] → out[0..4] (out starts empty, no prev accumulation) + lo = a0 * b0 + out[0] = lo + carry = unsignedMultiplyHigh(a0, b0) - val prev = out[i + j] - val s1 = prev + lo - val c1 = if (s1.toULong() < prev.toULong()) 1L else 0L - val s2 = s1 + carry - val c2 = if (s2.toULong() < s1.toULong()) 1L else 0L - out[i + j] = s2 - carry = hi + c1 + c2 - } - out[i + 4] = carry - } + lo = a0 * b1 + s = lo + carry + c1 = if (s.toULong() < lo.toULong()) 1L else 0L + out[1] = s + carry = unsignedMultiplyHigh(a0, b1) + c1 + + lo = a0 * b2 + s = lo + carry + c1 = if (s.toULong() < lo.toULong()) 1L else 0L + out[2] = s + carry = unsignedMultiplyHigh(a0, b2) + c1 + + lo = a0 * b3 + s = lo + carry + c1 = if (s.toULong() < lo.toULong()) 1L else 0L + out[3] = s + out[4] = unsignedMultiplyHigh(a0, b3) + c1 + + // Row 1: a1 × [b0,b1,b2,b3] accumulated into out[1..5] + lo = a1 * b0 + hi = unsignedMultiplyHigh(a1, b0) + prev = out[1] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + out[1] = s + carry = hi + c1 + + lo = a1 * b1 + hi = unsignedMultiplyHigh(a1, b1) + prev = out[2] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + s += carry + c2 = if (s.toULong() < carry.toULong()) 1L else 0L + out[2] = s + carry = hi + c1 + c2 + + lo = a1 * b2 + hi = unsignedMultiplyHigh(a1, b2) + prev = out[3] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + s += carry + c2 = if (s.toULong() < carry.toULong()) 1L else 0L + out[3] = s + carry = hi + c1 + c2 + + lo = a1 * b3 + hi = unsignedMultiplyHigh(a1, b3) + prev = out[4] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + s += carry + c2 = if (s.toULong() < carry.toULong()) 1L else 0L + out[4] = s + out[5] = hi + c1 + c2 + + // Row 2: a2 × [b0,b1,b2,b3] accumulated into out[2..6] + lo = a2 * b0 + hi = unsignedMultiplyHigh(a2, b0) + prev = out[2] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + out[2] = s + carry = hi + c1 + + lo = a2 * b1 + hi = unsignedMultiplyHigh(a2, b1) + prev = out[3] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + s += carry + c2 = if (s.toULong() < carry.toULong()) 1L else 0L + out[3] = s + carry = hi + c1 + c2 + + lo = a2 * b2 + hi = unsignedMultiplyHigh(a2, b2) + prev = out[4] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + s += carry + c2 = if (s.toULong() < carry.toULong()) 1L else 0L + out[4] = s + carry = hi + c1 + c2 + + lo = a2 * b3 + hi = unsignedMultiplyHigh(a2, b3) + prev = out[5] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + s += carry + c2 = if (s.toULong() < carry.toULong()) 1L else 0L + out[5] = s + out[6] = hi + c1 + c2 + + // Row 3: a3 × [b0,b1,b2,b3] accumulated into out[3..7] + lo = a3 * b0 + hi = unsignedMultiplyHigh(a3, b0) + prev = out[3] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + out[3] = s + carry = hi + c1 + + lo = a3 * b1 + hi = unsignedMultiplyHigh(a3, b1) + prev = out[4] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + s += carry + c2 = if (s.toULong() < carry.toULong()) 1L else 0L + out[4] = s + carry = hi + c1 + c2 + + lo = a3 * b2 + hi = unsignedMultiplyHigh(a3, b2) + prev = out[5] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + s += carry + c2 = if (s.toULong() < carry.toULong()) 1L else 0L + out[5] = s + carry = hi + c1 + c2 + + lo = a3 * b3 + hi = unsignedMultiplyHigh(a3, b3) + prev = out[6] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + s += carry + c2 = if (s.toULong() < carry.toULong()) 1L else 0L + out[6] = s + out[7] = hi + c1 + c2 } /** * Dedicated squaring: out = a² (512-bit result in LongArray(8)). * Exploits symmetry: 6 cross-products doubled + 4 diagonal = 10 multiplyHigh calls. + * Fully unrolled for ART JIT optimization. */ fun sqrWide( out: LongArray, a: LongArray, ) { - for (i in 0 until 8) out[i] = 0L + val a0 = a[0] + val a1 = a[1] + val a2 = a[2] + val a3 = a[3] + var lo: Long + var hi: Long + var prev: Long + var s: Long + var c1: Long + var c2: Long + var carry: Long + var v: Long - // Pass 1: cross-products a[i]*a[j] for i < j (single) - for (i in 0 until 4) { - var carry = 0L - val ai = a[i] - for (j in i + 1 until 4) { - val lo = ai * a[j] - val hi = unsignedMultiplyHigh(ai, a[j]) - val prev = out[i + j] - val s1 = prev + lo - val c1 = if (s1.toULong() < prev.toULong()) 1L else 0L - val s2 = s1 + carry - val c2 = if (s2.toULong() < s1.toULong()) 1L else 0L - out[i + j] = s2 - carry = hi + c1 + c2 - } - out[i + 4] = carry - } + // Pass 1: cross-products a[i]*a[j] for i < j (single, before doubling) + + // Row i=0: a0 × [a1, a2, a3] → out[1..4] + out[0] = 0L + lo = a0 * a1 + out[1] = lo + carry = unsignedMultiplyHigh(a0, a1) + + lo = a0 * a2 + s = lo + carry + c1 = if (s.toULong() < lo.toULong()) 1L else 0L + out[2] = s + carry = unsignedMultiplyHigh(a0, a2) + c1 + + lo = a0 * a3 + s = lo + carry + c1 = if (s.toULong() < lo.toULong()) 1L else 0L + out[3] = s + out[4] = unsignedMultiplyHigh(a0, a3) + c1 + + // Row i=1: a1 × [a2, a3] → accumulated into out[3..5] + lo = a1 * a2 + hi = unsignedMultiplyHigh(a1, a2) + prev = out[3] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + out[3] = s + carry = hi + c1 + + lo = a1 * a3 + hi = unsignedMultiplyHigh(a1, a3) + prev = out[4] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + s += carry + c2 = if (s.toULong() < carry.toULong()) 1L else 0L + out[4] = s + out[5] = hi + c1 + c2 + + // Row i=2: a2 × [a3] → accumulated into out[5..6] + lo = a2 * a3 + hi = unsignedMultiplyHigh(a2, a3) + prev = out[5] + s = prev + lo + c1 = if (s.toULong() < prev.toULong()) 1L else 0L + out[5] = s + out[6] = hi + c1 // Pass 2: double all cross-products (shift left by 1 bit) - var shiftCarry = 0L - for (i in 1 until 8) { - val v = out[i] - out[i] = (v shl 1) or shiftCarry - shiftCarry = v ushr 63 - } + v = out[1] + out[1] = v shl 1 + var shiftCarry = v ushr 63 + v = out[2] + out[2] = (v shl 1) or shiftCarry + shiftCarry = v ushr 63 + v = out[3] + out[3] = (v shl 1) or shiftCarry + shiftCarry = v ushr 63 + v = out[4] + out[4] = (v shl 1) or shiftCarry + shiftCarry = v ushr 63 + v = out[5] + out[5] = (v shl 1) or shiftCarry + shiftCarry = v ushr 63 + v = out[6] + out[6] = (v shl 1) or shiftCarry + shiftCarry = v ushr 63 + out[7] = shiftCarry // Pass 3: add diagonal products a[i]² - var dCarry = 0L - for (i in 0 until 4) { - val lo = a[i] * a[i] - val hi = unsignedMultiplyHigh(a[i], a[i]) - val pos = 2 * i - val s1 = out[pos] + lo - val c1 = if (s1.toULong() < out[pos].toULong()) 1L else 0L - val s2 = s1 + dCarry - val c2 = if (s2.toULong() < s1.toULong()) 1L else 0L - out[pos] = s2 + // i=0: a0², pos=0 + lo = a0 * a0 + hi = unsignedMultiplyHigh(a0, a0) + out[0] = lo // out[0] was 0 + s = out[1] + hi + c1 = if (s.toULong() < out[1].toULong()) 1L else 0L + out[1] = s + var dCarry = c1 - val s3 = out[pos + 1] + hi - val c3 = if (s3.toULong() < out[pos + 1].toULong()) 1L else 0L - val s4 = s3 + c1 + c2 - val c4 = if (s4.toULong() < s3.toULong()) 1L else 0L - out[pos + 1] = s4 - dCarry = c3 + c4 - } + // i=1: a1², pos=2 + lo = a1 * a1 + hi = unsignedMultiplyHigh(a1, a1) + s = out[2] + lo + c1 = if (s.toULong() < out[2].toULong()) 1L else 0L + s += dCarry + c2 = if (s.toULong() < dCarry.toULong()) 1L else 0L + out[2] = s + prev = out[3] + hi + val c3a = if (prev.toULong() < out[3].toULong()) 1L else 0L + prev += c1 + c2 + val c4a = if (prev.toULong() < (c1 + c2).toULong()) 1L else 0L + out[3] = prev + dCarry = c3a + c4a + + // i=2: a2², pos=4 + lo = a2 * a2 + hi = unsignedMultiplyHigh(a2, a2) + s = out[4] + lo + c1 = if (s.toULong() < out[4].toULong()) 1L else 0L + s += dCarry + c2 = if (s.toULong() < dCarry.toULong()) 1L else 0L + out[4] = s + prev = out[5] + hi + val c3b = if (prev.toULong() < out[5].toULong()) 1L else 0L + prev += c1 + c2 + val c4b = if (prev.toULong() < (c1 + c2).toULong()) 1L else 0L + out[5] = prev + dCarry = c3b + c4b + + // i=3: a3², pos=6 + lo = a3 * a3 + hi = unsignedMultiplyHigh(a3, a3) + s = out[6] + lo + c1 = if (s.toULong() < out[6].toULong()) 1L else 0L + s += dCarry + c2 = if (s.toULong() < dCarry.toULong()) 1L else 0L + out[6] = s + prev = out[7] + hi + val c3c = if (prev.toULong() < out[7].toULong()) 1L else 0L + prev += c1 + c2 + out[7] = prev } // ==================== Serialization ====================