perf: optimize secp256k1 field arithmetic and point operations for Android

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
This commit is contained in:
Claude
2026-04-07 23:14:11 +00:00
parent fc07090121
commit 1dd8edcd41
4 changed files with 586 additions and 174 deletions
@@ -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)
}
@@ -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.
@@ -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<AffinePoint>,
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). */
@@ -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 ====================