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 ====================