From adc9bbdf23611d3c8fe6595a43167e98ee97f9ff Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Tue, 21 Jul 2026 22:34:14 -0500 Subject: [PATCH 01/16] secp256k1: Correct comments in field 64 reduce. This corrects the referenced names in the proof in the comments of field64Reuce512 and tightens the bounds. While the bounds that were commented were sufficient to prove it is safe to discard the carry, using a tighter bound provides an even stronger proof. While here, also clarify in the field64Square512 comments that the impossibility of the p5 carry is formally proven in internal/proof to avoid potentially giving the wrong impression due omitting the case analysis is merely informal. --- dcrec/secp256k1/field64.go | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/dcrec/secp256k1/field64.go b/dcrec/secp256k1/field64.go index 9f4c58d3a..18651d92a 100644 --- a/dcrec/secp256k1/field64.go +++ b/dcrec/secp256k1/field64.go @@ -799,9 +799,10 @@ func field64Square512(t *[8]uint64, a *[4]uint64) { // The p5 carry is safe to discard because p5 + h13 + c ≤ 2^64 - 1 (where c // is the carry from p4 + l13). // - // A full proof involves case work that is omitted here, but the key point - // is that the only way the final add could have a carry is if all 3 of the - // following conditions were simultaneously true: + // A full proof involves case analysis that is omitted here since the + // impossibility of the carry is formally proven in internal/proofs, but the + // key point is that the only way the final add could have a carry is if all + // 3 of the following conditions were simultaneously true: // // 1) p5_old = 1 (the carry from the earlier chain, so ≤ 1) // 2) h13 = 2^64 - 2 (h13 ≤ 2^64 - 2 as proven previously) @@ -891,14 +892,15 @@ func field64Reduce512(r *[4]uint64, x *[8]uint64) { h, t0 = bits.Mul64(x[4], field64PrimeComplement) - // Note that since hi is the upper 64 bits of the product of two uint64s: - // h3 ≤ floor((2^64-1)^2 / 2^64) = 2^64 - 2 + // Note that since hi is the upper 64 bits of the product of a uint64 with + // c and c < 2^33: + // hi ≤ floor((2^64-1)(2^33 - 1) / 2^64) = 2^33 - 2 // - // Then, because c ≤ 1, a loose bound is: - // p4 ≤ h3 + 1 = 2^64 - 1 < 2^64 + // Then, because carry ≤ 1, a loose bound for h is: + // h ≤ hi + 1 = 2^33 - 1 < 2^64 // // Therefore, it is safe to discard the carry and the same applies to the - // next two limbs. + // next two limbs (second h and first t4). hi, lo = bits.Mul64(x[5], field64PrimeComplement) t1, carry = bits.Add64(lo, h, 0) h, _ = bits.Add64(hi, 0, carry) From 2379513a32fe738be36db44eeeaa920ac2ca883e Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Tue, 21 Jul 2026 22:34:19 -0500 Subject: [PATCH 02/16] secp256k1: Improve proof readability. This cleans up the Python Z3 proofs to express the propositions directly and remove a bunch of boilerplate. Currently, every proof creates a solver, asserts the logical negation of the proposition, and calls `check` . Likewise, each discarded carry proof repeats the same loop and solver setup. Also, proof failures are only printed, so the scripts continue executing instead of returning a non-zero exit status suitable for CI. This introduces `prove` and `prove_no_discarded_carries` helpers that encapsulate the common patterns and raise an assertion on failure. Rather than requiring each call site to express the negated proposition, `prove` accepts the proposition directly and internally shows that its negation is unsatisfiable. This hides the SMT-specific proof mechanism which allows the proof scripts to express the intended invariants more naturally. Finally, it updates the existing proof scripts to use the new helpers. --- .../proofs/field_4x64_reduce_prover.py | 20 +++---- .../internal/proofs/mul512_4x64_prover.py | 34 ++++-------- .../internal/proofs/square512_4x64_prover.py | 6 +-- .../internal/proofs/z3_proof_helpers.py | 52 +++++++++++++++---- 4 files changed, 59 insertions(+), 53 deletions(-) diff --git a/dcrec/secp256k1/internal/proofs/field_4x64_reduce_prover.py b/dcrec/secp256k1/internal/proofs/field_4x64_reduce_prover.py index d15241524..a5172e902 100644 --- a/dcrec/secp256k1/internal/proofs/field_4x64_reduce_prover.py +++ b/dcrec/secp256k1/internal/proofs/field_4x64_reduce_prover.py @@ -92,24 +92,18 @@ # ------- # Discarded carries are never set. -for idx, discarded in enumerate(discards): - s = Solver() - s.add(discarded != ZERO) - check(s, f"discarded carry {idx} != 0") +prove_no_discarded_carries(discards) # Top limb never exceeds the field complement after first reduction. -s = Solver() -s.add(UGT(t4_saved, field64PrimeComplement)) -check(s, "top limb after 1st reduction > complement") +# +# Since the complement is 33 bits, this proves the value does not exceed +# 256 + 33 = 289 bits after the first reduction. +prove(ULE(t4_saved, field64PrimeComplement), "top limb after 1st reduction > complement") # Value after second reduction is less than twice the prime (t < 2p). -s = Solver() t = Concat(t4, t3, t2, t1, t0) -s.add(UGE(t, twicePrime)) -check(s, "value after 2nd reduction >= 2p") +prove(ULT(t, twicePrime), "value after 2nd reduction >= 2p") # Fully reduced result is less than the prime (r < p). -s = Solver() r = Concat(r3, r2, r1, r0) -s.add(UGE(r, field64Prime)) -check(s, "fully reduced result >= prime") +prove(ULT(r, field64Prime), "fully reduced result >= prime") diff --git a/dcrec/secp256k1/internal/proofs/mul512_4x64_prover.py b/dcrec/secp256k1/internal/proofs/mul512_4x64_prover.py index 944864df0..13ceff260 100644 --- a/dcrec/secp256k1/internal/proofs/mul512_4x64_prover.py +++ b/dcrec/secp256k1/internal/proofs/mul512_4x64_prover.py @@ -99,27 +99,13 @@ # ------- # Discarded carries are never set. -for idx, discarded in enumerate(discards): - s = Solver() - s.add(discarded != ZERO) - check(s, f"discarded carry {idx} != 0") - -# Top limb after row 0 is max 2**64 - 2 (p4 ≤ 2**64 - 2). -s = Solver() -s.add(UGT(p4_saved_0, (1<<64) - 2)) -check(s, "top limb after row 0 > 2**64-2") - -# Limb 5 after row 1 is max 2**64 - 2 (q4 ≤ 2**64 - 2). -s = Solver() -s.add(UGT(q4_saved_1, (1<<64) - 2)) -check(s, "limb 5 after row 1 > 2**64-2") - -# Limb 5 after row 2 is max 2**64 - 2 (q4 ≤ 2**64 - 2). -s = Solver() -s.add(UGT(q4_saved_2, (1<<64) - 2)) -check(s, "limb 5 after row 2 > 2**64-2") - -# Limb 5 after row 3 is max 2**64 - 2 (q4 ≤ 2**64 - 2). -s = Solver() -s.add(UGT(q4_saved_3, (1<<64) - 2)) -check(s, "limb 5 after row 3 > 2**64-2") +prove_no_discarded_carries(discards) + +# Top limb of each of the four rows is max 2**64 - 2. +# +# For row 0, that corresponds to: (p4 ≤ 2**64 - 2) +# For rows 1-3, it corresponds to: (q4 ≤ 2**64 - 2) +prove(ULE(p4_saved_0, (1<<64) - 2), "top limb after row 0 > 2**64-2") +prove(ULE(q4_saved_1, (1<<64) - 2), "limb 5 after row 1 > 2**64-2") +prove(ULE(q4_saved_2, (1<<64) - 2), "limb 5 after row 2 > 2**64-2") +prove(ULE(q4_saved_3, (1<<64) - 2), "limb 5 after row 3 > 2**64-2") diff --git a/dcrec/secp256k1/internal/proofs/square512_4x64_prover.py b/dcrec/secp256k1/internal/proofs/square512_4x64_prover.py index cf4e01407..3248a93f8 100644 --- a/dcrec/secp256k1/internal/proofs/square512_4x64_prover.py +++ b/dcrec/secp256k1/internal/proofs/square512_4x64_prover.py @@ -13,7 +13,6 @@ a1 = BitVec('a1', 64) a2 = BitVec('a2', 64) a3 = BitVec('a3', 64) -a_full = Concat(a3, a2, a1, a0) # --------------- # Model the code. @@ -73,7 +72,4 @@ # ------- # Discarded carries are never set. -for idx, discarded in enumerate(discards): - s = Solver() - s.add(discarded != ZERO) - check(s, f"discarded carry {idx} != 0") +prove_no_discarded_carries(discards) diff --git a/dcrec/secp256k1/internal/proofs/z3_proof_helpers.py b/dcrec/secp256k1/internal/proofs/z3_proof_helpers.py index 02a5649c8..28ee409a2 100644 --- a/dcrec/secp256k1/internal/proofs/z3_proof_helpers.py +++ b/dcrec/secp256k1/internal/proofs/z3_proof_helpers.py @@ -30,7 +30,12 @@ def add64(x, y, cin): Return (sum, carry) exactly like bits.Add64. x,y : 64-bit bitvectors - cin : 0/1 as a 64-bit bitvector + cin : 0 or 1 as a 64-bit bitvector + + Note: + Although cin is a 64-bit bitvector, it should semantically represent a + single carry bit (0 or 1). Z3 handles arbitrary values correctly + (extracting only bit 0), but callers should maintain this invariant. """ assert x.size() == 64 assert y.size() == 64 @@ -45,7 +50,12 @@ def sub64(x, y, bin): Return (difference, borrow) exactly like bits.Sub64. x,y : 64-bit bitvectors - bin : 0/1 as a 64-bit bitvector + bin : 0 or 1 as a 64-bit bitvector + + Note: + Although bin is a 64-bit bitvector, it should semantically represent a + single borrow bit (0 or 1). Z3 handles arbitrary values correctly + (extracting only bit 0), but callers should maintain this invariant. """ assert x.size() == 64 assert y.size() == 64 @@ -57,14 +67,34 @@ def sub64(x, y, bin): return lo, borrow def check(s, reason): - check_result = s.check() - print(f"{reason}: {check_result}") + """Fail unless the proof succeeds.""" + result = s.check() + if result != unsat: + m = s.model() if result == sat else None + details = "" + if m is not None: + details = "\n " + ", ".join( + f"{d.name()}=0x{m[d].as_long():016x}" for d in m.decls() + ) + raise AssertionError(f"{reason}: {result}{details}") + print(f"{reason}: {result}") - if check_result == sat: - m = s.model() +def prove(proposition, reason, assumptions=()): + """ + Prove the proposition is valid by showing its negation is unsatisfiable. - for v in [x0,x1,x2,x3,x4,x5,x6,x7]: - if m[v] != None: - print(f"{v} = 0x{m[v].as_long():016x}") - else: - print(f"{v} = {m[v]}") + The assumptions encode additional assumptions proven elsehwere. + """ + s = Solver() + for a in assumptions: + s.add(a) + s.add(Not(proposition)) + check(s, reason) + +def prove_no_discarded_carries(discards, label="discarded carry"): + """ + Prove every carry in discards is always zero. This is to say, discarding it + never loses information. + """ + for idx, discarded in enumerate(discards): + prove(discarded == ZERO, f"{label} {idx} != 0") From 383f7a041f601c718ac2c941b8c2f3d232fc2c6e Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Tue, 21 Jul 2026 22:34:20 -0500 Subject: [PATCH 03/16] secp256k1: Use QF_BV for Z3 proofs. The proofs operate entirely within the quantifier-free bitvector (QF_BV) logic, so there is no need to use the generic solver. This updates `prove` to construct a QF_BV solver that is specialized for the logic being proved directly which avoids the overhead of selecting from more general solving strategies. --- dcrec/secp256k1/internal/proofs/z3_proof_helpers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dcrec/secp256k1/internal/proofs/z3_proof_helpers.py b/dcrec/secp256k1/internal/proofs/z3_proof_helpers.py index 28ee409a2..089d75267 100644 --- a/dcrec/secp256k1/internal/proofs/z3_proof_helpers.py +++ b/dcrec/secp256k1/internal/proofs/z3_proof_helpers.py @@ -85,7 +85,7 @@ def prove(proposition, reason, assumptions=()): The assumptions encode additional assumptions proven elsehwere. """ - s = Solver() + s = SolverFor("QF_BV") for a in assumptions: s.add(a) s.add(Not(proposition)) From 250306d29c6a9787d359b211869f0d2257fe2209 Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Tue, 21 Jul 2026 22:34:21 -0500 Subject: [PATCH 04/16] secp256k1: Compactify Z3 proofs. This makes use of BitVecVals to specify multiple vectors at the same time which both improves readability and makes the proofs slightly more compact. It also makes the field64 reducer proof more closely match the source code which accesses the limb of x by array index. --- .../proofs/field_4x64_reduce_prover.py | 47 ++++++++++--------- .../internal/proofs/mul512_4x64_prover.py | 10 +--- .../internal/proofs/square512_4x64_prover.py | 5 +- 3 files changed, 27 insertions(+), 35 deletions(-) diff --git a/dcrec/secp256k1/internal/proofs/field_4x64_reduce_prover.py b/dcrec/secp256k1/internal/proofs/field_4x64_reduce_prover.py index a5172e902..c0d7ef582 100644 --- a/dcrec/secp256k1/internal/proofs/field_4x64_reduce_prover.py +++ b/dcrec/secp256k1/internal/proofs/field_4x64_reduce_prover.py @@ -9,21 +9,20 @@ # Inputs. # ------- -x0 = BitVec('x0', 64) -x1 = BitVec('x1', 64) -x2 = BitVec('x2', 64) -x3 = BitVec('x3', 64) -x4 = BitVec('x4', 64) -x5 = BitVec('x5', 64) -x6 = BitVec('x6', 64) -x7 = BitVec('x7', 64) +x = BitVecs('x0 x1 x2 x3 x4 x5 x6 x7', 64) +# NOTE: The code in this section reuses the variable names from the Go code +# being proven for easy cross referencing. However, the comments refer to them +# by shorter mathematical variables for clarity. + +# Define P. field64Prime0 = BitVecVal(0xfffffffefffffc2f, 64) field64Prime1 = BitVecVal(0xffffffffffffffff, 64) field64Prime2 = BitVecVal(0xffffffffffffffff, 64) field64Prime3 = BitVecVal(0xffffffffffffffff, 64) field64Prime = Concat(field64Prime3, field64Prime2, field64Prime1, field64Prime0) +# Define 2P. twicePrime0 = BitVecVal(0xfffffffdfffff85e, 64) twicePrime1 = BitVecVal(0xffffffffffffffff, 64) twicePrime2 = BitVecVal(0xffffffffffffffff, 64) @@ -31,6 +30,7 @@ twicePrime4 = ONE twicePrime = Concat(twicePrime4, twicePrime3, twicePrime2, twicePrime1, twicePrime0) +# Define c = 2**256 - P. field64PrimeComplement = BitVecVal(2**32+977, 64) # --------------- @@ -40,33 +40,34 @@ discards = [] # first reduction: 512 bits -> 289 bits. -h, t0 = mul64(x4, field64PrimeComplement) +h, t0 = mul64(x[4], field64PrimeComplement) -hi, lo = mul64(x5, field64PrimeComplement) +hi, lo = mul64(x[5], field64PrimeComplement) t1, carry = add64(lo, h, ZERO) h, discarded = add64(hi, ZERO, carry) discards.append(discarded) -hi, lo = mul64(x6, field64PrimeComplement) +hi, lo = mul64(x[6], field64PrimeComplement) t2, carry = add64(lo, h, ZERO) h, discarded = add64(hi, ZERO, carry) discards.append(discarded) -hi, lo = mul64(x7, field64PrimeComplement) +hi, lo = mul64(x[7], field64PrimeComplement) t3, carry = add64(lo, h, ZERO) t4, discarded = add64(hi, ZERO, carry) discards.append(discarded) -t0, carry = add64(t0, x0, ZERO) -t1, carry = add64(t1, x1, carry) -t2, carry = add64(t2, x2, carry) -t3, carry = add64(t3, x3, carry) -t4 += carry +t0, carry = add64(t0, x[0], ZERO) +t1, carry = add64(t1, x[1], carry) +t2, carry = add64(t2, x[2], carry) +t3, carry = add64(t3, x[3], carry) +t4, discarded = add64(t4, ZERO, carry) +discards.append(discarded) # Save for analysis. t4_saved = t4 -# second reduction: 289 bits -> t < 2p +# second reduction: 289 bits -> t < 2P h, t4 = mul64(t4, field64PrimeComplement) t0, carry = add64(t0, t4, ZERO) @@ -76,7 +77,7 @@ t4 = carry -# final reduction: t < 2p -> t < p +# final reduction: t < 2P -> t < P s0, borrow = sub64(t0, field64Prime0, ZERO) s1, borrow = sub64(t1, field64Prime1, borrow) s2, borrow = sub64(t2, field64Prime2, borrow) @@ -100,10 +101,10 @@ # 256 + 33 = 289 bits after the first reduction. prove(ULE(t4_saved, field64PrimeComplement), "top limb after 1st reduction > complement") -# Value after second reduction is less than twice the prime (t < 2p). +# Value after second reduction is less than twice the prime (t < 2P). t = Concat(t4, t3, t2, t1, t0) -prove(ULT(t, twicePrime), "value after 2nd reduction >= 2p") +prove(ULT(t, twicePrime), "value after 2nd reduction >= 2P") -# Fully reduced result is less than the prime (r < p). +# Fully reduced result is less than the prime (r < P). r = Concat(r3, r2, r1, r0) -prove(ULT(r, field64Prime), "fully reduced result >= prime") +prove(ULT(r, field64Prime), "fully reduced result >= P") diff --git a/dcrec/secp256k1/internal/proofs/mul512_4x64_prover.py b/dcrec/secp256k1/internal/proofs/mul512_4x64_prover.py index 13ceff260..955f174ca 100644 --- a/dcrec/secp256k1/internal/proofs/mul512_4x64_prover.py +++ b/dcrec/secp256k1/internal/proofs/mul512_4x64_prover.py @@ -9,14 +9,8 @@ # Inputs. # ------- -a0 = BitVec('a0', 64) -a1 = BitVec('a1', 64) -a2 = BitVec('a2', 64) -a3 = BitVec('a3', 64) -b0 = BitVec('b0', 64) -b1 = BitVec('b1', 64) -b2 = BitVec('b2', 64) -b3 = BitVec('b3', 64) +a0, a1, a2, a3 = BitVecs('a0 a1 a2 a3', 64) +b0, b1, b2, b3 = BitVecs('b0 b1 b2 b3', 64) # --------------- # Model the code. diff --git a/dcrec/secp256k1/internal/proofs/square512_4x64_prover.py b/dcrec/secp256k1/internal/proofs/square512_4x64_prover.py index 3248a93f8..0cb5da8e9 100644 --- a/dcrec/secp256k1/internal/proofs/square512_4x64_prover.py +++ b/dcrec/secp256k1/internal/proofs/square512_4x64_prover.py @@ -9,10 +9,7 @@ # Inputs. # ------- -a0 = BitVec('a0', 64) -a1 = BitVec('a1', 64) -a2 = BitVec('a2', 64) -a3 = BitVec('a3', 64) +a0, a1, a2, a3 = BitVecs('a0 a1 a2 a3', 64) # --------------- # Model the code. From 2b1ec4c7f26c5cbedcf9c47ccee03ed472e5ac99 Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Tue, 21 Jul 2026 22:34:23 -0500 Subject: [PATCH 05/16] secp256k1: Add formal functional field redux proof. The existing formal proof for the 4x64 field reduction establishes the intermediate bounds and verifies that no carries are discarded. This extends it to include a full end-to-end functional correctness proof that shows the implementation computes the correct reduction modulo the field prime for all 512-bit inputs. --- .../proofs/field_4x64_reduce_prover.py | 235 ++++++++++++++++++ .../internal/proofs/z3_proof_helpers.py | 5 + 2 files changed, 240 insertions(+) diff --git a/dcrec/secp256k1/internal/proofs/field_4x64_reduce_prover.py b/dcrec/secp256k1/internal/proofs/field_4x64_reduce_prover.py index c0d7ef582..2c6a7aeb4 100644 --- a/dcrec/secp256k1/internal/proofs/field_4x64_reduce_prover.py +++ b/dcrec/secp256k1/internal/proofs/field_4x64_reduce_prover.py @@ -33,6 +33,29 @@ # Define c = 2**256 - P. field64PrimeComplement = BitVecVal(2**32+977, 64) +# -------------------------------------------------------------------------- +# Prove self-consistency of the constants before trusting anything that uses +# them. +# -------------------------------------------------------------------------- + +# Prove twicePrime = 2P. +field64Prime320 = ZeroExt(64, field64Prime) +prove(twicePrime == field64Prime320 << 1, "twicePrime != 2P") + +# -------------------------------------------------------------------------- +# Lemma R-identity: P + c == 2**256 +# +# This lemma justifies replacing multiples of 2**256 with multiples of c when +# working modulo P: +# +# P + c ≡ 2**256 (mod P) +# => 0 + c ≡ 2**256 (mod P) +# => c ≡ 2**256 (mod P) +# -------------------------------------------------------------------------- +field64PrimeComplement320 = ZeroExt(256, field64PrimeComplement) +prove(field64Prime320+field64PrimeComplement320 == BitVecVal(1, 320)<<256, + "Lemma 1: P+c != 2**256") + # --------------- # Model the code. # --------------- @@ -108,3 +131,215 @@ # Fully reduced result is less than the prime (r < P). r = Concat(r3, r2, r1, r0) prove(ULT(r, field64Prime), "fully reduced result >= P") + +# ----------------------------------------------------------------------------- +# Prove functional congruence: r ≡ x (mod P). +# +# The checks above establish bounds and the safety of discarded carries, +# self-consistency of the constants, and the reduction identity (lemma +# R-identity). +# +# The following lemmas establish the arithmetic identity of each transformation +# stage. Together with the previously proven bounds and reduction identity, +# they imply the overall congruence r ≡ x (mod P). +# +# A direct assertion such as t == x_lo + (x >> 256)*P is intractable because +# it introduces a monolithic 256x256 multiplier that appears nowhere in the +# modeled code, so the SAT solver is forced into bit-level equivalence +# checking between two structurally unrelated multiplier circuits. +# +# Instead, because every accumulation block in the implementation is a fixed +# carry-propagation network over the *outputs* of bits.Mul64, whether or not +# those 128-bit inputs are actually products is irrelevant to the network. Each +# block satisfies its per-block identity for ALL 128-bit values. So each block +# identity is proven here as a universal lemma over unconstrained variables, +# which contains no multipliers at all and reduces to adder-network equivalence +# that the solver dispatches quickly. +# +# The Go function, field64Reduce512, computes r in three stages: +# +# 1) First reduction: Fold the high 256 bits of x (x4..x7) into the low 256 bits +# (x0..x3). This produces t_1 = x_lo + x_hi*c which fits in 289 bits (5 +# 64-bit limbs with the 5th limb small). +# 2) Second reduction: Fold t_1's 5th limb back in the same way. This produces +# t_2 = t_1_lo + t_1_hi*c, where t_2 < 2P. +# 3) Final reduction: One conditional subtraction of P which is valid because +# t_2 < 2p. +# +# The lemmas below chain together to prove r ≡ t2 ≡ t1 ≡ x (mod P) and 0 ≤ r < +# P which is exactly the definition of x ≡ r (mod P). +def prove_functional_congruence_lemmas(): + """ + Prove r == x (mod P) for field64Reduce512. + + x is the 512-bit input (8 limbs x0..x7) and r is the fully-reduced 256-bit + output (4 limbs). + """ + + # Use a wide enough congruence width that nothing in a 512-bit input ever + # wraps around. This is comfortably wider than the largest quantity + # of x_hi*c < 2**289. + CONGRUENCE_WIDTH = 512 + + def cwidth(v): + """Zero-extend v to CONGRUENCE_WIDTH bits.""" + return ZeroExt(CONGRUENCE_WIDTH - v.size(), v) + + def concat(*limbs): + """ + Concatenate the provided big-endian 64-bit limbs into a single + CONGRUENCE_WIDTH-bit value: + + limb0*2**(64*num_limbs) + limb1*2**(64*(num_limbs-1)) + ... + """ + return cwidth(Concat(*limbs)) + + def climbs(limbs): + """ + Concatenate the provided iterable little-endian 64-bit limbs into a + single CONGRUENCE_WIDTH-bit value: + + limbs[0] + limbs[1]*2**64 + limbs[2]*2**128 + ... + """ + return cwidth(Concat(*reversed(limbs))) + + # Universal lemma variable shared by R-open, R-interior, and R-second. + # + # prod models whatever 128-bit product bits.Mul64(_, c) happens to return at + # that place in the real code. It is NOT tied to any specific x[j] or t4. + # hi/lo are its two halves. + prod = BitVec("prod", 128) + hi, lo = split128(prod) + + # ------------------------------------------------------------------------- + # Lemma R-open: The opening row of the first reduction (the x4 line). + # + # h, t0 = bits.Mul64(x[4], c) + # + # No addition happens on this line at all. It just splits the product. + # + # This lemma is really just confirming that packing lo and hi back together + # recreates the original product, which sounds trivial, but it matters + # because it's the base case that R-interior's telescoping argument (below) + # builds on. + # + # It establishes that after processing x4, the running total (h, t0) exactly + # equals x4*c with nothing lost. + # ------------------------------------------------------------------------- + prove(concat(hi,lo) == cwidth(prod), "R-open: (h,t0) != x4*c") + + # ------------------------------------------------------------------------- + # Lemma R-interior: The shape shared by the x5, x6, and x7 rows. + # + # The x7 row simply names its final output t4 instead of h. + # + # hi, lo = bits.Mul64(x[j], c) + # t_j, carry = bits.Add64(lo, h, 0) + # h, _ = bits.Add64(hi, 0, carry) + # + # h_in models the incoming running-total limb carried forward from the + # previous row (R-open's h, or a previous R-interior's h_out). + # + # One carry is discarded and assumed zero. This fact is proven above by the + # discard checks for the real x5/x6/x7 rows, so this does not smuggle in a + # new assumption. + # ------------------------------------------------------------------------- + h_in = BitVec("h_in", 64) + t_j, carry = add64(lo, h_in, ZERO) + h_out, discarded = add64(hi, ZERO, carry) + prove( + concat(h_out, t_j) == cwidth(h_in) + cwidth(prod), + "R-interior: (t_j,h_out) != h_in + xj*c", + assumptions=[discarded == ZERO], + ) + + # ------------------------------------------------------------------------- + # Lemma R-tail: Fold x0..x3 into running total after the four rows above. + # + # t0, carry = bits.Add64(t0, x[0], 0) + # t1, carry = bits.Add64(t1, x[1], carry) + # t2, carry = bits.Add64(t2, x[2], carry) + # t3, carry = bits.Add64(t3, x[3], carry) + # t4 += carry + # + # b0..b4 model the incoming running total in t0..t4. + # + # The final line is modeled here as add64 with its own carry-out assumed + # zero. This fact is proven above by the discard checks for the real total, + # so this does not smuggle in a new assumption. + # ------------------------------------------------------------------------- + b = BitVecs("b0 b1 b2 b3 b4", 64) + v = [None]*5 + v[0], carry = add64(b[0], x[0], ZERO) + v[1], carry = add64(b[1], x[1], carry) + v[2], carry = add64(b[2], x[2], carry) + v[3], carry = add64(b[3], x[3], carry) + v[4], discarded = add64(b[4], ZERO, carry) + prove( + climbs(v) == climbs(b) + climbs(x[:4]), + "R-tail: post-fold total != state + x_lo", + assumptions=[discarded == ZERO], + ) + + # ------------------------------------------------------------------------- + # Lemma R-second: The second reduction. + # + # h, t4 = bits.Mul64(t4, c) + # t0, carry = bits.Add64(t0, t4, 0) + # t1, carry = bits.Add64(t1, h, carry) + # t2, carry = bits.Add64(t2, 0, carry) + # t3, carry = bits.Add64(t3, 0, carry) + # t4 = carry + # + # prod is reused here for t4_old*c since it's still just a free 128-bit + # variable and nothing ties it to being x4*c specifically. + # + # d0..d3 model t0..t3 going into this stage. + # + # No carries are discarded and the final carry becomes the new t4, so no + # assumption is needed. + # ------------------------------------------------------------------------- + d = BitVecs("d0 d1 d2 d3", 64) + w = [None] * 5 + w[0], carry = add64(d[0], lo, ZERO) + w[1], carry = add64(d[1], hi, carry) + w[2], carry = add64(d[2], ZERO, carry) + w[3], carry = add64(d[3], ZERO, carry) + w[4] = carry + prove( + climbs(w) == climbs(d) + cwidth(prod), + "R-second: 2nd reduction != state + t4*c", + ) + + # ----------------------------------------------------------------- + # Lemma R-final: the conditional-subtraction final reduction. + # + # s0, borrow = bits.Sub64(t0, p0, 0) + # ... + # _, borrow = bits.Sub64(t4, 0, borrow) + # r[i] = constantTimeSelect64(borrow, t[i], s[i]) + # + # z0..z4 model t0..t4 going into this stage. + # + # This is only valid because t < 2P as otherwise t - P could still exceed P + # and the 4-limb result wouldn't represent it. This fact is proven above + # for the real total, so this does not smuggle in a new assumption. + # ----------------------------------------------------------------- + p = cwidth(field64Prime) + twop = cwidth(twicePrime) + z = BitVecs("z0 z1 z2 z3 z4", 64) + s = [None]*4 + s[0], borrow = sub64(z[0], field64Prime0, ZERO) + s[1], borrow = sub64(z[1], field64Prime1, borrow) + s[2], borrow = sub64(z[2], field64Prime2, borrow) + s[3], borrow = sub64(z[3], field64Prime3, borrow) + _, borrow = sub64(z[4], ZERO, borrow) + t = climbs(z) + r = If(borrow == ONE, t, climbs(s)) + prove( + r == If(ULT(t, p), t, t - p), + "R-final: conditional-subtract result != t mod P", + assumptions=[ULT(t, twop)], + ) + +prove_functional_congruence_lemmas() diff --git a/dcrec/secp256k1/internal/proofs/z3_proof_helpers.py b/dcrec/secp256k1/internal/proofs/z3_proof_helpers.py index 089d75267..db10679db 100644 --- a/dcrec/secp256k1/internal/proofs/z3_proof_helpers.py +++ b/dcrec/secp256k1/internal/proofs/z3_proof_helpers.py @@ -66,6 +66,11 @@ def sub64(x, y, bin): ONE, ZERO) return lo, borrow +def split128(product): + """Return (high, low) 64-bit halves from a 128-bit value.""" + assert product.size() == 128 + return Extract(127, 64, product), Extract(63, 0, product) + def check(s, reason): """Fail unless the proof succeeds.""" result = s.check() From fb39615ca0edd83af6fe239c6a7a98dafb71b70a Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Thu, 9 Jul 2026 01:28:21 -0500 Subject: [PATCH 06/16] secp256k1: Add caveats to documentation. This updates the README.md and doc.go files to include some important caveats about the constant time and memory clearing behavior. Specifically, while the behavior has been manually verified by carefully reviewing the generated assembly as of the most recent version of Go, it is impossible to guarantee the compiler will not change in the future. --- dcrec/secp256k1/README.md | 19 +++++++++++++++++-- dcrec/secp256k1/doc.go | 15 +++++++++++++++ 2 files changed, 32 insertions(+), 2 deletions(-) diff --git a/dcrec/secp256k1/README.md b/dcrec/secp256k1/README.md index 53f383911..4fcfa1861 100644 --- a/dcrec/secp256k1/README.md +++ b/dcrec/secp256k1/README.md @@ -16,8 +16,8 @@ structures and functions for working with public and private secp256k1 keys. In addition, subpackages are provided to produce, verify, parse, and serialize ECDSA signatures and EC-Schnorr-DCRv0 (a custom Schnorr-based signature scheme -specific to Decred) signatures. See the README.md files in the relevant sub -packages for more details about those aspects. +specific to Decred) signatures. See the README.md files in the relevant +subpackages for more details about those aspects. ## Features @@ -95,6 +95,21 @@ secp256k1 elliptic curve cryptography. The secp256k1 domain parameters are specified in [SEC 2: Recommended Elliptic Curve Domain Parameters](https://www.secg.org/sec2-v2.pdf). +## Caveats + +Constant time implementations are designed to eliminate data dependent timing +variation at the algorithmic level. This property depends on the Go compiler +continuing to compile the relevant operations into constant-time machine +instructions. + +This package also consistently uses temporary stack allocated buffers which it +attempts to clear in order to reduce the likelihood of sensitive data lingering +in memory, but Go does not guarantee such writes are never optimized away. + +Significant effort has been made to ensure correctness including careful review +of the generated assembly, but it is impossible to guarantee exact behavior of +future versions of the Go compiler. + ## secp256k1 use in Decred At the time of this writing, the primary public key cryptography in widespread diff --git a/dcrec/secp256k1/doc.go b/dcrec/secp256k1/doc.go index d946650d7..69121d116 100644 --- a/dcrec/secp256k1/doc.go +++ b/dcrec/secp256k1/doc.go @@ -50,6 +50,21 @@ use optimized secp256k1 elliptic curve cryptography. Finally, a comprehensive suite of tests is provided to provide a high level of quality assurance. +## Caveats + +Constant time implementations are designed to eliminate data dependent timing +variation at the algorithmic level. This property depends on the Go compiler +continuing to compile the relevant operations into constant time machine +instructions. + +This package also consistently uses temporary stack allocated buffers which it +attempts to clear in order to reduce the likelihood of sensitive data lingering +in memory, but Go does not guarantee such writes are never optimized away. + +Significant effort has been made to ensure correctness, including careful review +of the generated assembly, but it is impossible to guarantee exact behavior of +future versions of the Go compiler. + # Use of secp256k1 in Decred At the time of this writing, the primary public key cryptography in widespread From 779e55b2328f213bce4cde08d01f3de3341504fb Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Thu, 9 Jul 2026 03:35:47 -0500 Subject: [PATCH 07/16] secp256k1: Add scalar negation algebraic tests. This adds a new test helper for some algebraic properties of scalar negation and calls it from the test that handles specific edge cases as well as the randomized test. --- dcrec/secp256k1/modnscalar_test.go | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/dcrec/secp256k1/modnscalar_test.go b/dcrec/secp256k1/modnscalar_test.go index 251708a6f..b894e9b4c 100644 --- a/dcrec/secp256k1/modnscalar_test.go +++ b/dcrec/secp256k1/modnscalar_test.go @@ -979,6 +979,26 @@ func TestModNScalarSquareRandom(t *testing.T) { } } +// checkModNScalarNegateProps checks algebraic properties of scalar negation. +func checkModNScalarNegateProps(t *testing.T, val *ModNScalar) { + t.Helper() + + negation := new(ModNScalar).NegateVal(val) + + // Ensure involution produces the original value. That is -(-x) == x. + involution := new(ModNScalar).Set(negation).Negate() + if !involution.Equals(val) { + t.Errorf("mismatched involution -- got: %v, want: %v", involution, + val) + } + + // Ensure x + (-x) == 0. + sum := new(ModNScalar).Add2(val, negation) + if !sum.IsZero() { + t.Errorf("x + (-x) != 0 -- got: %v, want: 0", sum) + } +} + // TestModNScalarNegate ensures that negating scalars works as expected for edge // cases. func TestModNScalarNegate(t *testing.T) { @@ -1051,6 +1071,9 @@ func TestModNScalarNegate(t *testing.T) { result2, expected) continue } + + // Ensure algebraic properties with negation hold. + checkModNScalarNegateProps(t, s) } } @@ -1086,6 +1109,9 @@ func TestModNScalarNegateRandom(t *testing.T) { "big int result: %x\nscalar result %v", bigIntVal, modNVal, bigIntResult, modNValResult) } + + // Ensure algebraic properties with negation hold. + assertModNScalarNegateProperties(t, modNVal) } } From c4c734c98e1164148af9d4af103197f2b1707eb2 Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Thu, 9 Jul 2026 04:10:05 -0500 Subject: [PATCH 08/16] secp256k1: Add scalar inverse algebraic tests. This adds a new test helper for some algebraic properties of modular multiplicative inverses of scalars and calls it from the test that handles specific edge cases as well as the randomized test. --- dcrec/secp256k1/modnscalar_test.go | 34 +++++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/dcrec/secp256k1/modnscalar_test.go b/dcrec/secp256k1/modnscalar_test.go index b894e9b4c..e4a9307b6 100644 --- a/dcrec/secp256k1/modnscalar_test.go +++ b/dcrec/secp256k1/modnscalar_test.go @@ -1111,7 +1111,33 @@ func TestModNScalarNegateRandom(t *testing.T) { } // Ensure algebraic properties with negation hold. - assertModNScalarNegateProperties(t, modNVal) + checkModNScalarNegateProps(t, modNVal) + } +} + +// checkModNScalarInverseProps checks algebraic properties of the +// multiplicative modular inverse of scalars. +func checkModNScalarInverseProps(t *testing.T, val *ModNScalar) { + t.Helper() + + // Exception for 0 which does not have a modular multiplicative inverse. + if val.IsZero() { + return + } + + inverse := new(ModNScalar).InverseValNonConst(val) + + // Ensure x * x^-1 ≡ 1 (mod N). + product := new(ModNScalar).Mul2(val, inverse) + one := new(ModNScalar).SetInt(1) + if !product.Equals(one) { + t.Errorf("x * x^-1 != 1 (mod N) -- got: %v, want: %v", product, one) + } + + // Ensure (x^-1)^-1 == x. + doubleInv := new(ModNScalar).InverseValNonConst(inverse) + if !doubleInv.Equals(val) { + t.Errorf("(x^-1)^-1 != x (mod N) -- got: %v, want: %v", doubleInv, val) } } @@ -1189,6 +1215,9 @@ func TestModNScalarInverseNonConst(t *testing.T) { result2, expected) continue } + + // Ensure algebraic properties with inverses hold. + checkModNScalarInverseProps(t, s) } } @@ -1224,6 +1253,9 @@ func TestModNScalarInverseNonConstRandom(t *testing.T) { "big int result: %x\nscalar result %v", bigIntVal, modNVal, bigIntResult, modNValResult) } + + // Ensure algebraic properties with inverses hold. + checkModNScalarInverseProps(t, modNVal) } } From 2a9a517e7e5e84741e05e77cee009fd987ab07b7 Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Thu, 9 Jul 2026 04:30:46 -0500 Subject: [PATCH 09/16] secp256k1: Add scalar half order algebraic tests. This adds a new test helper for some algebraic properties of the scalar half order check and calls it from the test that handles specific edge cases as well as the randomized test. --- dcrec/secp256k1/modnscalar_test.go | 30 +++++++++++++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/dcrec/secp256k1/modnscalar_test.go b/dcrec/secp256k1/modnscalar_test.go index e4a9307b6..448378e0c 100644 --- a/dcrec/secp256k1/modnscalar_test.go +++ b/dcrec/secp256k1/modnscalar_test.go @@ -1259,6 +1259,27 @@ func TestModNScalarInverseNonConstRandom(t *testing.T) { } } +// checkModNScalarHalfOrderProps checks algebraic properties of the half order +// comparison. +func checkModNScalarHalfOrderProps(t *testing.T, val *ModNScalar) { + t.Helper() + + // Exception for 0 which is the only scalar where negation does not produce + // a distinct value. + if val.IsZero() { + return + } + + // Ensure negating a scalar flips whether it is over half order. + // + // x > N/2 implies N-x < N/2 and vice versa. + negation := new(ModNScalar).NegateVal(val) + if val.IsOverHalfOrder() == negation.IsOverHalfOrder() { + t.Errorf("negation did not flip half order -- val: %v, negation: %v", + val, negation) + } +} + // TestModNScalarIsOverHalfOrder ensures that scalars report whether or not they // exceeed the half order works as expected for edge cases. func TestModNScalarIsOverHalfOrder(t *testing.T) { @@ -1305,12 +1326,16 @@ func TestModNScalarIsOverHalfOrder(t *testing.T) { }} for _, test := range tests { - result := mustModNScalarWithOverflow(test.in).IsOverHalfOrder() + s := mustModNScalarWithOverflow(test.in) + result := s.IsOverHalfOrder() if result != test.expected { t.Errorf("%s: unexpected result -- got: %v, want: %v", test.name, result, test.expected) continue } + + // Ensure algebraic properties with the half order hold. + checkModNScalarHalfOrderProps(t, s) } } @@ -1344,5 +1369,8 @@ func TestModNScalarIsOverHalfOrderRandom(t *testing.T) { "in: %v\nbig int result: %v\nscalar result %v", bigIntVal, modNVal, bigIntResult, modNValResult) } + + // Ensure algebraic properties with the half order hold. + checkModNScalarHalfOrderProps(t, modNVal) } } From ab27bafec6aa8baca1d6b9343e7af55b33bc25fe Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Thu, 9 Jul 2026 01:28:29 -0500 Subject: [PATCH 10/16] secp256k1: Implement 4x64 ModNScalar. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This significantly optimizes the ModNScalar implementation by rewriting it to use saturated 4x64 arithmetic similar to the new 4x64 field arithmetic. The new implementation retains all of the important properties such as being constant time and only using temporary buffers on the stack that are cleared. Of particular interest is the new modular reduction implementation. It moves away from the 96-bit accumulate design and instead manually tracks and discards carries inline which provides a major speed advantage. Due to the complexity and potential mistakes in the arithmetic, all of the intermediate bounds and carry assumptions documented that the code relies upon have been formally verified to be correct. The formal proof is in a separate commit. The implementation is primarily targeted at 64-bit since it is by far the most common architecture in modern use, however, as the included benchmarks below show, it is faster for almost all operations on 32-bit architectures as well. The few operations where it is a bit slower on 32-bit architectures are not used anywhere near as frequently as the ones that get 10-30% better performance, so the end result is the overall implementation is notably faster on 32 bit too. For 64-bit architectures, every operation is significantly faster. For example, all of the primary operations that get heavy use get ~73-94% better performance. Finally, the tests have been updated to use values specifically crafted to exercise the new limbs versus the old values that were crafted for the old limbs. Comparison for 32-bit platforms: $ benchstat beforescalar_x86.txt afterscalar_x86.txt name old time/op new time/op delta ----------------------------------------------------------------------------------- ModNScalar 22.9ns ± 4% 23.2ns ± 5% ~ (p=0.184 n=10+10) ModNScalarZero 1.36ns ± 3% 6.59ns ± 6% +385.95% (p=0.000 n=9+9) ModNScalarIsZero 0.32ns ± 4% 0.34ns ±14% ~ (p=0.143 n=10+10) ModNScalarEquals 3.69ns ± 4% 0.30ns ± 7% -91.78% (p=0.000 n=10+10) ModNScalarAdd 25.1ns ± 6% 22.5ns ± 5% -10.40% (p=0.000 n=9+10) ModNScalarMul 257ns ± 5% 223ns ± 5% -13.07% (p=0.000 n=10+10) ModNScalarSquare 267ns ± 5% 172ns ± 6% -35.56% (p=0.000 n=10+10) ModNScalarNegate 12.8ns ± 9% 14.5ns ± 2% +13.39% (p=0.000 n=10+10) ModNScalarInverse 686ns ± 7% 623ns ± 3% -9.24% (p=0.000 n=10+9) ModNScalarIsOverHalfOrder 8.92ns ± 6% 6.55ns ± 3% -26.59% (p=0.000 n=10+10) Comparison for 64-bit platforms: $ benchstat beforescalar_x64.txt afterscalar_x64.txt name old time/op new time/op delta ModNScalar 12.1ns ± 1% 3.5ns ± 4% -71.25% (p=0.000 n=7+10) ModNScalarZero 0.64ns ± 9% 0.61ns ± 5% -4.08% (p=0.043 n=10+10) ModNScalarIsZero 0.31ns ± 8% 0.31ns ± 6% ~ (p=0.000 n=10+10) ModNScalarEquals 2.63ns ± 5% 0.33ns ± 9% -87.47% (p=0.000 n=10+10) ModNScalarAdd 15.1ns ± 4% 4.1ns ± 4% -73.25% (p=0.000 n=10+10) ModNScalarMul 214ns ± 6% 29ns ± 6% -86.30% (p=0.000 n=10+10) ModNScalarSquare 215ns ± 4% 23ns ± 6% -89.39% (p=0.000 n=10+10) ModNScalarNegate 5.78ns ± 4% 2.99ns ±10% -48.33% (p=0.000 n=10+10) ModNScalarInverse 452ns ± 4% 408ns ± 6% -9.89% (p=0.000 n=10+10) ModNScalarIsOverHalfOrder 5.52ns ± 6% 0.29ns ± 5% -94.66% (p=0.000 n=10+9) --- dcrec/secp256k1/common.go | 7 +- dcrec/secp256k1/curve.go | 54 +- dcrec/secp256k1/curve_test.go | 20 +- dcrec/secp256k1/field_test.go | 4 +- dcrec/secp256k1/modnscalar.go | 1366 ++++++++++++---------------- dcrec/secp256k1/modnscalar_test.go | 473 +++++----- 6 files changed, 881 insertions(+), 1043 deletions(-) diff --git a/dcrec/secp256k1/common.go b/dcrec/secp256k1/common.go index b45e013bd..9ff782bf4 100644 --- a/dcrec/secp256k1/common.go +++ b/dcrec/secp256k1/common.go @@ -15,9 +15,10 @@ func constantTimeEq64(a, b uint64) uint32 { return uint32(((t | -t) >> 63) ^ 1) } -// constantTimeNotEq returns 1 if a != b or 0 otherwise in constant time. -func constantTimeNotEq(a, b uint32) uint32 { - return ^uint32((uint64(a^b)-1)>>63) & 1 +// constantTimeNotEq64 returns 1 if a != b or 0 otherwise in constant time. +func constantTimeNotEq64(a, b uint64) uint32 { + t := a ^ b + return uint32(((t | -t) >> 63)) } // constantTimeLess returns 1 if a < b or 0 otherwise in constant time. diff --git a/dcrec/secp256k1/curve.go b/dcrec/secp256k1/curve.go index 84c6d0c41..750b57a32 100644 --- a/dcrec/secp256k1/curve.go +++ b/dcrec/secp256k1/curve.go @@ -720,16 +720,6 @@ func mulAdd64Carry(digit1, digit2, m, c uint64) (hi, lo uint64) { // because it is used for replacing division with multiplication and thus the // intermediate results must be done via a field extension to a larger field. func mul512Rsh320Round(n1, n2 *ModNScalar) ModNScalar { - // Convert n1 and n2 to base 2^64 digits. - n1Digit0 := uint64(n1.n[0]) | uint64(n1.n[1])<<32 - n1Digit1 := uint64(n1.n[2]) | uint64(n1.n[3])<<32 - n1Digit2 := uint64(n1.n[4]) | uint64(n1.n[5])<<32 - n1Digit3 := uint64(n1.n[6]) | uint64(n1.n[7])<<32 - n2Digit0 := uint64(n2.n[0]) | uint64(n2.n[1])<<32 - n2Digit1 := uint64(n2.n[2]) | uint64(n2.n[3])<<32 - n2Digit2 := uint64(n2.n[4]) | uint64(n2.n[5])<<32 - n2Digit3 := uint64(n2.n[6]) | uint64(n2.n[7])<<32 - // Compute the full 512-bit product n1*n2. var r0, r1, r2, r3, r4, r5, r6, r7, c uint64 @@ -738,40 +728,40 @@ func mul512Rsh320Round(n1, n2 *ModNScalar) ModNScalar { // // Note that r0 is ignored because it is not needed to compute the higher // terms and it is shifted out below anyway. - c, _ = bits.Mul64(n2Digit0, n1Digit0) - c, r1 = mulAdd64(n2Digit0, n1Digit1, c) - c, r2 = mulAdd64(n2Digit0, n1Digit2, c) - r4, r3 = mulAdd64(n2Digit0, n1Digit3, c) + c, _ = bits.Mul64(n2.n[0], n1.n[0]) + c, r1 = mulAdd64(n2.n[0], n1.n[1], c) + c, r2 = mulAdd64(n2.n[0], n1.n[2], c) + r4, r3 = mulAdd64(n2.n[0], n1.n[3], c) // Terms resulting from the product of the second digit of the second number // by all digits of the first number. // // Note that r1 is ignored because it is no longer needed to compute the // higher terms and it is shifted out below anyway. - c, _ = mulAdd64(n2Digit1, n1Digit0, r1) - c, r2 = mulAdd64Carry(n2Digit1, n1Digit1, r2, c) - c, r3 = mulAdd64Carry(n2Digit1, n1Digit2, r3, c) - r5, r4 = mulAdd64Carry(n2Digit1, n1Digit3, r4, c) + c, _ = mulAdd64(n2.n[1], n1.n[0], r1) + c, r2 = mulAdd64Carry(n2.n[1], n1.n[1], r2, c) + c, r3 = mulAdd64Carry(n2.n[1], n1.n[2], r3, c) + r5, r4 = mulAdd64Carry(n2.n[1], n1.n[3], r4, c) // Terms resulting from the product of the third digit of the second number // by all digits of the first number. // // Note that r2 is ignored because it is no longer needed to compute the // higher terms and it is shifted out below anyway. - c, _ = mulAdd64(n2Digit2, n1Digit0, r2) - c, r3 = mulAdd64Carry(n2Digit2, n1Digit1, r3, c) - c, r4 = mulAdd64Carry(n2Digit2, n1Digit2, r4, c) - r6, r5 = mulAdd64Carry(n2Digit2, n1Digit3, r5, c) + c, _ = mulAdd64(n2.n[2], n1.n[0], r2) + c, r3 = mulAdd64Carry(n2.n[2], n1.n[1], r3, c) + c, r4 = mulAdd64Carry(n2.n[2], n1.n[2], r4, c) + r6, r5 = mulAdd64Carry(n2.n[2], n1.n[3], r5, c) // Terms resulting from the product of the fourth digit of the second number // by all digits of the first number. // // Note that r3 is ignored because it is no longer needed to compute the // higher terms and it is shifted out below anyway. - c, _ = mulAdd64(n2Digit3, n1Digit0, r3) - c, r4 = mulAdd64Carry(n2Digit3, n1Digit1, r4, c) - c, r5 = mulAdd64Carry(n2Digit3, n1Digit2, r5, c) - r7, r6 = mulAdd64Carry(n2Digit3, n1Digit3, r6, c) + c, _ = mulAdd64(n2.n[3], n1.n[0], r3) + c, r4 = mulAdd64Carry(n2.n[3], n1.n[1], r4, c) + c, r5 = mulAdd64Carry(n2.n[3], n1.n[2], r5, c) + r7, r6 = mulAdd64Carry(n2.n[3], n1.n[3], r6, c) // At this point the upper 256 bits of the full 512-bit product n1*n2 are in // r4..r7 (recall the low order results were discarded as noted above). @@ -796,14 +786,10 @@ func mul512Rsh320Round(n1, n2 *ModNScalar) ModNScalar { // less than the group order given the group order is > 2^255 and the // maximum possible value of the result is 2^192. var result ModNScalar - result.n[0] = uint32(r0) - result.n[1] = uint32(r0 >> 32) - result.n[2] = uint32(r1) - result.n[3] = uint32(r1 >> 32) - result.n[4] = uint32(r2) - result.n[5] = uint32(r2 >> 32) - result.n[6] = uint32(r3) - result.n[7] = uint32(r3 >> 32) + result.n[0] = r0 + result.n[1] = r1 + result.n[2] = r2 + result.n[3] = r3 return result } diff --git a/dcrec/secp256k1/curve_test.go b/dcrec/secp256k1/curve_test.go index b1b4bec00..d3eb497f2 100644 --- a/dcrec/secp256k1/curve_test.go +++ b/dcrec/secp256k1/curve_test.go @@ -877,28 +877,16 @@ func TestScalarBaseMultJacobian(t *testing.T) { // modNBitLen returns the minimum number of bits required to represent the mod n // scalar. The result is 0 when the value is 0. func modNBitLen(s *ModNScalar) uint16 { - if w := s.n[7]; w > 0 { - return uint16(bits.Len32(w)) + 224 - } - if w := s.n[6]; w > 0 { - return uint16(bits.Len32(w)) + 192 - } - if w := s.n[5]; w > 0 { - return uint16(bits.Len32(w)) + 160 - } - if w := s.n[4]; w > 0 { - return uint16(bits.Len32(w)) + 128 - } if w := s.n[3]; w > 0 { - return uint16(bits.Len32(w)) + 96 + return uint16(bits.Len64(w)) + 192 } if w := s.n[2]; w > 0 { - return uint16(bits.Len32(w)) + 64 + return uint16(bits.Len64(w)) + 128 } if w := s.n[1]; w > 0 { - return uint16(bits.Len32(w)) + 32 + return uint16(bits.Len64(w)) + 64 } - return uint16(bits.Len32(s.n[0])) + return uint16(bits.Len64(s.n[0])) } // checkLambdaDecomposition returns an error if the provided decomposed scalars diff --git a/dcrec/secp256k1/field_test.go b/dcrec/secp256k1/field_test.go index a1200b8dd..41a5a6175 100644 --- a/dcrec/secp256k1/field_test.go +++ b/dcrec/secp256k1/field_test.go @@ -1307,7 +1307,7 @@ func TestFieldMulInt(t *testing.T) { in2: 5, expected: "6cf4eb20f2447c77657fccb172d38c0aa91ea4ac446dc641fa463a6b5091fba7", }, { - name: "random sampling #3", + name: "random sampling #4", in1: "fb4529be3e027a3d1587d8a500b72f2d312e3577340ef5175f96d113be4c2ceb", in2: 8, expected: "da294df1f013d1e8ac3ec52805b979698971abb9a077a8bafcb688a4f261820f", @@ -1369,7 +1369,7 @@ func TestFieldMul(t *testing.T) { in2: "3", expected: "0", }, { - name: "secp256k1 prime * 3", + name: "secp256k1 prime * 8", in1: "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2e", in2: "8", expected: "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc27", diff --git a/dcrec/secp256k1/modnscalar.go b/dcrec/secp256k1/modnscalar.go index 9405a987f..6ccaf6306 100644 --- a/dcrec/secp256k1/modnscalar.go +++ b/dcrec/secp256k1/modnscalar.go @@ -5,8 +5,10 @@ package secp256k1 import ( + "encoding/binary" "encoding/hex" "math/big" + "math/bits" "sync" ) @@ -18,81 +20,76 @@ import ( // https://cacr.uwaterloo.ca/hac/ // Many elliptic curve operations require working with scalars in a finite field -// characterized by the order of the group underlying the secp256k1 curve. -// Given this precision is larger than the biggest available native type, -// obviously some form of bignum math is needed. This code implements -// specialized fixed-precision field arithmetic rather than relying on an -// arbitrary-precision arithmetic package such as math/big for dealing with the -// math modulo the group order since the size is known. As a result, rather -// large performance gains are achieved by taking advantage of many -// optimizations not available to arbitrary-precision arithmetic and generic -// modular arithmetic algorithms. +// defined by the order of the group underlying the secp256k1 curve. Since the +// group order exceeds the width of the largest available native type, some form +// of multi-precision arithmetic is required. This code implements specialized +// fixed-precision field arithmetic rather than relying on an +// arbitrary-precision arithmetic package such as math/big for performing +// arithmetic modulo the fixed group order. As a result, significant +// performance gains are achieved by taking advantage of many optimizations not +// available to arbitrary-precision arithmetic and generic modular arithmetic +// algorithms. // -// There are various ways to internally represent each element. For example, -// the most obvious representation would be to use an array of 4 uint64s (64 -// bits * 4 = 256 bits). However, that representation suffers from the fact -// that there is no native Go type large enough to handle the intermediate -// results while adding or multiplying two 64-bit numbers. +// There are various ways to internally represent each element with their own +// pros and cons. Some common representations are: // -// Given the above, this implementation represents the field elements as 8 -// uint32s with each word (array entry) treated as base 2^32. This was chosen -// because most systems at the current time are 64-bit (or at least have 64-bit -// registers available for specialized purposes such as MMX) so the intermediate -// results can typically be done using a native register (and using uint64s to -// avoid the need for additional half-word arithmetic) - +// - 8 uint32s with base 2^32 limbs (aka 8x32: 32 bits * 8 = 256 bits) +// - 10 uint32s with base 2^26 limbs (aka 10x26: 26 bits * 10 = 260 bits) +// - 4 uint64s with base 2^64 limbs (aka 4x64: 64 bits * 4 = 256 bits) +// - 5 uint64s with base 2^52 limbs (aka 5x52: 52 bits * 5 = 260 bits) +// +// The primary tradeoff is complexity versus performance and the optimal choice +// also largely depends on the available hardware capabilities. +// +// For example, 5x52 and 10x26 allow performing several operations in a row +// before carry propagation or modular reduction becomes necessary since there +// are additional bits available in each limb, but that comes at the cost of +// manual magnitude tracking, multiple non-canonical representations of each +// number, and periodic normalization to propagate the accumulated carries and +// perform the modular reduction all at once. +// +// Conversely, 8x32 and 4x64 involve less complexity, are simpler to use, have +// fewer limbs to manage, and only have a single canonical representation for +// each number, but that comes at the cost of requiring carry propagation and +// modular reduction for every operation and larger intermediate results. +// +// The requirement for larger intermediate results is particularly notable for +// uint64s because there is no native Go type large enough to handle the 128-bit +// intermediate results while adding and multiplying the 64-bit limbs. +// +// This implementation historically used 8x32 for that reason, coupled with the +// fact that Go did not reliably make use of hardware-supported wide arithmetic, +// to ensure all intermediate results fit cleanly into uint64s. +// +// However, most systems are now 64-bit and almost all have instructions capable +// of handling the 128-bit intermediate results without the need for additional +// half-word arithmetic. Further, modern versions of Go automatically use those +// instructions on hardware that supports it. +// +// Given the above, this implementation represents the field elements as 4 +// uint64s with each limb (array entry) treated as base 2^64 (aka 4x64) and +// keeps the value reduced at all times. const ( - // These fields provide convenient access to each of the words of the + // These fields provide convenient access to each of the limbs of the // secp256k1 curve group order N to improve code readability. // // The group order of the curve per [SECG] is: - // 0xffffffff ffffffff ffffffff fffffffe baaedce6 af48a03b bfd25e8c d0364141 - // - // nolint: dupword - orderWordZero uint32 = 0xd0364141 - orderWordOne uint32 = 0xbfd25e8c - orderWordTwo uint32 = 0xaf48a03b - orderWordThree uint32 = 0xbaaedce6 - orderWordFour uint32 = 0xfffffffe - orderWordFive uint32 = 0xffffffff - orderWordSix uint32 = 0xffffffff - orderWordSeven uint32 = 0xffffffff - - // These fields provide convenient access to each of the words of the two's + // 0xfffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141 + orderLimb0 uint64 = 0xbfd25e8cd0364141 + orderLimb1 uint64 = 0xbaaedce6af48a03b + orderLimb2 uint64 = 0xfffffffffffffffe + orderLimb3 uint64 = 0xffffffffffffffff + + // These fields provide convenient access to each of the limbs of the two's // complement of the secp256k1 curve group order N to improve code // readability. // // The two's complement of the group order is: // 0x00000000 00000000 00000000 00000001 45512319 50b75fc4 402da173 2fc9bebf - orderComplementWordZero uint32 = (^orderWordZero) + 1 - orderComplementWordOne uint32 = ^orderWordOne - orderComplementWordTwo uint32 = ^orderWordTwo - orderComplementWordThree uint32 = ^orderWordThree - // orderComplementWordFour uint32 = ^orderWordFour // unused - // orderComplementWordFive uint32 = ^orderWordFive // unused - // orderComplementWordSix uint32 = ^orderWordSix // unused - // orderComplementWordSeven uint32 = ^orderWordSeven // unused - - // These fields provide convenient access to each of the words of the - // secp256k1 curve group order N / 2 to improve code readability and avoid - // the need to recalculate them. - // - // The half order of the secp256k1 curve group is: - // 0x7fffffff ffffffff ffffffff ffffffff 5d576e73 57a4501d dfe92f46 681b20a0 - // - // nolint: dupword - halfOrderWordZero uint32 = 0x681b20a0 - halfOrderWordOne uint32 = 0xdfe92f46 - halfOrderWordTwo uint32 = 0x57a4501d - halfOrderWordThree uint32 = 0x5d576e73 - halfOrderWordFour uint32 = 0xffffffff - halfOrderWordFive uint32 = 0xffffffff - halfOrderWordSix uint32 = 0xffffffff - halfOrderWordSeven uint32 = 0x7fffffff - - // uint32Mask is simply a mask with all bits set for a uint32 and is used to - // improve the readability of the code. - uint32Mask = 0xffffffff + orderComplementLimb0 uint64 = (^orderLimb0) + 1 + orderComplementLimb1 uint64 = ^orderLimb1 + orderComplementLimb2 uint64 = ^orderLimb2 + // orderComplementLimb3 uint64 = ^orderLimb3 // unused ) var ( @@ -112,39 +109,37 @@ var ( // around if absolutely needed. For example, subtraction can be performed by // adding the negation. // -// Should it be absolutely necessary, conversion to the standard library -// math/big.Int can be accomplished by using the Bytes method, slicing the -// resulting fixed-size array, and feeding it to big.Int.SetBytes. However, -// that should typically be avoided when possible as conversion to big.Ints -// requires allocations, is not constant time, and is slower when working modulo -// the group order. +// Should it be absolutely necessary, conversion to a standard library [big.Int] +// can be accomplished by using [ModNScalar.Bytes], slicing the resulting +// fixed-size array, and feeding it to [big.Int.SetBytes]. However, that should +// typically be avoided when possible as conversion to a [big.Int] requires +// allocations, is not constant time, and is slower when working modulo the +// group order. type ModNScalar struct { - // The scalar is represented as 8 32-bit integers in base 2^32. + // The scalar is represented as 4 64-bit integers in base 2^64. // // The following depicts the internal representation: - // --------------------------------------------------------- - // | n[7] | n[6] | ... | n[0] | - // | 32 bits | 32 bits | ... | 32 bits | - // | Mult: 2^(32*7) | Mult: 2^(32*6) | ... | Mult: 2^(32*0) | - // --------------------------------------------------------- + // -------------------------------------------------------------------- + // | n[3] | n[2] | n[1] | n[0] | + // | 64 bits | 64 bits | 64 bits | 64 bits | + // | Mult: 2^(64*3) | Mult: 2^(64*2) | Mult: 2^(64*1) | Mult: 2^(64*0) | + // -------------------------------------------------------------------- // - // For example, consider the number 2^87 + 2^42 + 1. It would be + // For example, consider the number 2^135 + 2^87 + 1. It would be // represented as: // n[0] = 1 - // n[1] = 2^10 - // n[2] = 2^23 - // n[3..7] = 0 - // - // The full 256-bit value is then calculated by looping i from 7..0 and - // doing sum(n[i] * 2^(32i)) like so: - // n[7] * 2^(32*7) = 0 * 2^224 = 0 - // n[6] * 2^(32*6) = 0 * 2^192 = 0 - // ... - // n[2] * 2^(32*2) = 2^23 * 2^64 = 2^87 - // n[1] * 2^(32*1) = 2^10 * 2^32 = 2^42 - // n[0] * 2^(32*0) = 1 * 2^0 = 1 - // Sum: 0 + 0 + ... + 2^87 + 2^42 + 1 = 2^87 + 2^42 + 1 - n [8]uint32 + // n[1] = 2^23 + // n[2] = 2^7 + // n[3] = 0 + // + // The full 256-bit value is then calculated by looping i from 3..0 and + // doing sum(n[i] * 2^(64i)) like so: + // n[3] * 2^(64*3) = 0 * 2^192 = 0 + // n[2] * 2^(64*2) = 2^7 * 2^128 = 2^135 + // n[1] * 2^(64*1) = 2^23 * 2^64 = 2^87 + // n[0] * 2^(64*0) = 1 * 2^0 = 1 + // Sum: 0 + 2^135 + 2^87 + 1 = 2^135 + 2^87 + 1 + n [4]uint64 } // String returns the scalar as a human-readable hex string. @@ -169,14 +164,7 @@ func (s *ModNScalar) Set(val *ModNScalar) *ModNScalar { // already set to zero. This function can be useful to clear an existing scalar // for reuse. func (s *ModNScalar) Zero() { - s.n[0] = 0 - s.n[1] = 0 - s.n[2] = 0 - s.n[3] = 0 - s.n[4] = 0 - s.n[5] = 0 - s.n[6] = 0 - s.n[7] = 0 + s.n = [4]uint64{} } // IsZeroBit returns 1 when the scalar is equal to zero or 0 otherwise in @@ -187,16 +175,13 @@ func (s *ModNScalar) Zero() { // operations require a numeric value. See [ModNScalar.IsZero] for the version // that returns a bool. func (s *ModNScalar) IsZeroBit() uint32 { - // The scalar can only be zero if no bits are set in any of the words. - bits := s.n[0] | s.n[1] | s.n[2] | s.n[3] | s.n[4] | s.n[5] | s.n[6] | s.n[7] - return constantTimeEq(bits, 0) + return constantTimeEq64(s.n[0]|s.n[1]|s.n[2]|s.n[3], 0) } // IsZero returns whether or not the scalar is equal to zero in constant time. func (s *ModNScalar) IsZero() bool { // The scalar can only be zero if no bits are set in any of the words. - bits := s.n[0] | s.n[1] | s.n[2] | s.n[3] | s.n[4] | s.n[5] | s.n[6] | s.n[7] - return bits == 0 + return (s.n[0] | s.n[1] | s.n[2] | s.n[3]) == 0 } // SetInt sets the scalar to the passed integer in constant time. This is a @@ -206,79 +191,10 @@ func (s *ModNScalar) IsZero() bool { // The scalar is returned to support chaining. This enables syntax like: // s := new(ModNScalar).SetInt(2).Mul(s2) so that s = 2 * s2. func (s *ModNScalar) SetInt(ui uint32) *ModNScalar { - s.Zero() - s.n[0] = ui + s.n = [4]uint64{uint64(ui), 0, 0, 0} return s } -// overflows determines if the current scalar is greater than or equal to the -// group order in constant time and returns 1 if it is or 0 otherwise. -func (s *ModNScalar) overflows() uint32 { - // The intuition here is that the scalar is greater than the group order if - // one of the higher individual words is greater than corresponding word of - // the group order and all higher words in the scalar are equal to their - // corresponding word of the group order. Since this type is modulo the - // group order, being equal is also an overflow back to 0. - // - // Note that the words 5, 6, and 7 are all the max uint32 value, so there is - // no need to test if those individual words of the scalar exceeds them, - // hence, only equality is checked for them. - highWordsEqual := constantTimeEq(s.n[7], orderWordSeven) - highWordsEqual &= constantTimeEq(s.n[6], orderWordSix) - highWordsEqual &= constantTimeEq(s.n[5], orderWordFive) - overflow := highWordsEqual & constantTimeGreater(s.n[4], orderWordFour) - highWordsEqual &= constantTimeEq(s.n[4], orderWordFour) - overflow |= highWordsEqual & constantTimeGreater(s.n[3], orderWordThree) - highWordsEqual &= constantTimeEq(s.n[3], orderWordThree) - overflow |= highWordsEqual & constantTimeGreater(s.n[2], orderWordTwo) - highWordsEqual &= constantTimeEq(s.n[2], orderWordTwo) - overflow |= highWordsEqual & constantTimeGreater(s.n[1], orderWordOne) - highWordsEqual &= constantTimeEq(s.n[1], orderWordOne) - overflow |= highWordsEqual & constantTimeGreaterOrEq(s.n[0], orderWordZero) - - return overflow -} - -// reduce256 reduces the current scalar modulo the group order in accordance -// with the overflows parameter in constant time. The overflows parameter -// specifies whether or not the scalar is known to be greater than the group -// order and MUST either be 1 in the case it is or 0 in the case it is not for a -// correct result. -func (s *ModNScalar) reduce256(overflows uint32) { - // Notice that since s < 2^256 < 2N (where N is the group order), the max - // possible number of reductions required is one. Therefore, in the case a - // reduction is needed, it can be performed with a single subtraction of N. - // Also, recall that subtraction is equivalent to addition by the two's - // complement while ignoring the carry. - // - // When s >= N, the overflows parameter will be 1. Conversely, it will be 0 - // when s < N. Thus multiplying by the overflows parameter will either - // result in 0 or the multiplicand itself. - // - // Combining the above along with the fact that s + 0 = s, the following is - // a constant time implementation that works by either adding 0 or the two's - // complement of N as needed. - // - // The final result will be in the range 0 <= s < N as expected. - overflows64 := uint64(overflows) - c := uint64(s.n[0]) + overflows64*uint64(orderComplementWordZero) - s.n[0] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(s.n[1]) + overflows64*uint64(orderComplementWordOne) - s.n[1] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(s.n[2]) + overflows64*uint64(orderComplementWordTwo) - s.n[2] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(s.n[3]) + overflows64*uint64(orderComplementWordThree) - s.n[3] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(s.n[4]) + overflows64 // * 1 - s.n[4] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(s.n[5]) // + overflows64 * 0 - s.n[5] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(s.n[6]) // + overflows64 * 0 - s.n[6] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(s.n[7]) // + overflows64 * 0 - s.n[7] = uint32(c & uint32Mask) -} - // SetBytes interprets the provided array as a 256-bit big-endian unsigned // integer, reduces it modulo the group order, sets the scalar to the result, // and returns either 1 if it was reduced (aka it overflowed) or 0 otherwise in @@ -288,23 +204,41 @@ func (s *ModNScalar) reduce256(overflows uint32) { // from a bool to numeric value in constant time and many constant-time // operations require a numeric value. func (s *ModNScalar) SetBytes(b *[32]byte) uint32 { - // Pack the 256 total bits across the 8 uint32 words. This could be done + // Pack the 256 total bits across the 4 uint64 words. This could be done // with a for loop, but benchmarks show this unrolled version is about 2 // times faster than the variant that uses a loop. - s.n[0] = uint32(b[31]) | uint32(b[30])<<8 | uint32(b[29])<<16 | uint32(b[28])<<24 - s.n[1] = uint32(b[27]) | uint32(b[26])<<8 | uint32(b[25])<<16 | uint32(b[24])<<24 - s.n[2] = uint32(b[23]) | uint32(b[22])<<8 | uint32(b[21])<<16 | uint32(b[20])<<24 - s.n[3] = uint32(b[19]) | uint32(b[18])<<8 | uint32(b[17])<<16 | uint32(b[16])<<24 - s.n[4] = uint32(b[15]) | uint32(b[14])<<8 | uint32(b[13])<<16 | uint32(b[12])<<24 - s.n[5] = uint32(b[11]) | uint32(b[10])<<8 | uint32(b[9])<<16 | uint32(b[8])<<24 - s.n[6] = uint32(b[7]) | uint32(b[6])<<8 | uint32(b[5])<<16 | uint32(b[4])<<24 - s.n[7] = uint32(b[3]) | uint32(b[2])<<8 | uint32(b[1])<<16 | uint32(b[0])<<24 - - // The value might be >= N, so reduce it as required and return whether or - // not it was reduced. - needsReduce := s.overflows() - s.reduce256(needsReduce) - return needsReduce + s.n[0] = binary.BigEndian.Uint64(b[24:32]) + s.n[1] = binary.BigEndian.Uint64(b[16:24]) + s.n[2] = binary.BigEndian.Uint64(b[8:16]) + s.n[3] = binary.BigEndian.Uint64(b[0:8]) + + // Since s < 2^256 < 2N (where N is the secp256k1 group order), the max + // possible number of reductions required is one. Therefore, in the case a + // reduction is needed, it can be performed with a single subtraction of N. + // + // Since N must only conditionally be subtracted when s ≥ N, the following + // handles it in constant time by always calculating t = s - N and selecting + // the correct case via a constant time select. + + // Subtract N with borrow propagation. borrow is set iff s < N. + // + // In other words, the input overflowed (≥ N) when s - N does NOT borrow. + // + // t = s - N + var t0, t1, t2, t3, borrow uint64 + t0, borrow = bits.Sub64(s.n[0], orderLimb0, 0) + t1, borrow = bits.Sub64(s.n[1], orderLimb1, borrow) + t2, borrow = bits.Sub64(s.n[2], orderLimb2, borrow) + t3, borrow = bits.Sub64(s.n[3], orderLimb3, borrow) + + // Constant-time select. + // + // Set s = s when s < N (aka borrow is set). Otherwise s = t = s - N. + s.n[0] = constantTimeSelect64(borrow, s.n[0], t0) + s.n[1] = constantTimeSelect64(borrow, s.n[1], t1) + s.n[2] = constantTimeSelect64(borrow, s.n[2], t2) + s.n[3] = constantTimeSelect64(borrow, s.n[3], t3) + return uint32(1 - borrow) } // zeroArray32 zeroes the provided 32-byte buffer. @@ -324,6 +258,8 @@ func zeroArray32(b *[32]byte) { // size or it is acceptable to use this function with the described truncation // and overflow behavior. func (s *ModNScalar) SetByteSlice(b []byte) bool { + // Always copy a total of 32 bytes regardless of the input length to avoid + // introducing data-dependent timing. var b32 [32]byte b = b[:constantTimeMin(uint32(len(b)), 32)] copy(b32[:], b32[:32-len(b)]) @@ -337,49 +273,21 @@ func (s *ModNScalar) SetByteSlice(b []byte) bool { // into the passed byte slice in constant time. The target slice must have at // least 32 bytes available or it will panic. // -// There is a similar function, PutBytes, which unpacks the scalar into a -// 32-byte array directly. This version is provided since it can be useful to -// write directly into part of a larger buffer without needing a separate -// allocation. +// There is a similar function, [ModNScalar.PutBytes], which unpacks the scalar +// into a 32-byte array directly. This version is provided since it can be +// useful to write directly into part of a larger buffer without needing a +// separate allocation. // // Preconditions: // - The target slice MUST have at least 32 bytes available func (s *ModNScalar) PutBytesUnchecked(b []byte) { - // Unpack the 256 total bits from the 8 uint32 words. This could be done - // with a for loop, but benchmarks show this unrolled version is about 2 + // Unpack the 256 total bits from the 4 uint64 words. This could be done + // with a for loop, but benchmarks show this unrolled version is about 3 // times faster than the variant which uses a loop. - b[31] = byte(s.n[0]) - b[30] = byte(s.n[0] >> 8) - b[29] = byte(s.n[0] >> 16) - b[28] = byte(s.n[0] >> 24) - b[27] = byte(s.n[1]) - b[26] = byte(s.n[1] >> 8) - b[25] = byte(s.n[1] >> 16) - b[24] = byte(s.n[1] >> 24) - b[23] = byte(s.n[2]) - b[22] = byte(s.n[2] >> 8) - b[21] = byte(s.n[2] >> 16) - b[20] = byte(s.n[2] >> 24) - b[19] = byte(s.n[3]) - b[18] = byte(s.n[3] >> 8) - b[17] = byte(s.n[3] >> 16) - b[16] = byte(s.n[3] >> 24) - b[15] = byte(s.n[4]) - b[14] = byte(s.n[4] >> 8) - b[13] = byte(s.n[4] >> 16) - b[12] = byte(s.n[4] >> 24) - b[11] = byte(s.n[5]) - b[10] = byte(s.n[5] >> 8) - b[9] = byte(s.n[5] >> 16) - b[8] = byte(s.n[5] >> 24) - b[7] = byte(s.n[6]) - b[6] = byte(s.n[6] >> 8) - b[5] = byte(s.n[6] >> 16) - b[4] = byte(s.n[6] >> 24) - b[3] = byte(s.n[7]) - b[2] = byte(s.n[7] >> 8) - b[1] = byte(s.n[7] >> 16) - b[0] = byte(s.n[7] >> 24) + binary.BigEndian.PutUint64(b[0:8], s.n[3]) + binary.BigEndian.PutUint64(b[8:16], s.n[2]) + binary.BigEndian.PutUint64(b[16:24], s.n[1]) + binary.BigEndian.PutUint64(b[24:32], s.n[0]) } // PutBytes unpacks the scalar to a 32-byte big-endian value using the passed @@ -419,11 +327,8 @@ func (s *ModNScalar) IsOdd() bool { func (s *ModNScalar) Equals(val *ModNScalar) bool { // Xor only sets bits when they are different, so the two scalars can only // be the same if no bits are set after xoring each word. - bits := (s.n[0] ^ val.n[0]) | (s.n[1] ^ val.n[1]) | (s.n[2] ^ val.n[2]) | - (s.n[3] ^ val.n[3]) | (s.n[4] ^ val.n[4]) | (s.n[5] ^ val.n[5]) | - (s.n[6] ^ val.n[6]) | (s.n[7] ^ val.n[7]) - - return bits == 0 + return ((s.n[0] ^ val.n[0]) | (s.n[1] ^ val.n[1]) | (s.n[2] ^ val.n[2]) | + (s.n[3] ^ val.n[3])) == 0 } // Add2 adds the passed two scalars together modulo the group order in constant @@ -431,27 +336,42 @@ func (s *ModNScalar) Equals(val *ModNScalar) bool { // // The scalar is returned to support chaining. This enables syntax like: // s3.Add2(s, s2).AddInt(1) so that s3 = s + s2 + 1. -func (s *ModNScalar) Add2(val1, val2 *ModNScalar) *ModNScalar { - c := uint64(val1.n[0]) + uint64(val2.n[0]) - s.n[0] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(val1.n[1]) + uint64(val2.n[1]) - s.n[1] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(val1.n[2]) + uint64(val2.n[2]) - s.n[2] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(val1.n[3]) + uint64(val2.n[3]) - s.n[3] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(val1.n[4]) + uint64(val2.n[4]) - s.n[4] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(val1.n[5]) + uint64(val2.n[5]) - s.n[5] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(val1.n[6]) + uint64(val2.n[6]) - s.n[6] = uint32(c & uint32Mask) - c = (c >> 32) + uint64(val1.n[7]) + uint64(val2.n[7]) - s.n[7] = uint32(c & uint32Mask) - - // The result is now 256 bits, but it might still be >= N, so use the - // existing normal reduce method for 256-bit values. - s.reduce256(uint32(c>>32) + s.overflows()) +func (s *ModNScalar) Add2(a, b *ModNScalar) *ModNScalar { + // Since both values are already in the range 0 ≤ val < N (where N is the + // secp256k1 group order), the maximum possible result is < 2N - 1. So a + // maximum of one subtraction of N is required in the worst case. + // + // Since N must only conditionally be subtracted when a+b ≥ N, the following + // handles it in constant time by calculating both t = a+b and u = a+b - N + // and selecting the correct case via a constant time select. + + // Add with carry propagation. overflow is set iff t = a+b ≥ 2^256. + // + // t = a + b + var t0, t1, t2, t3, overflow, carry uint64 + t0, carry = bits.Add64(a.n[0], b.n[0], 0) + t1, carry = bits.Add64(a.n[1], b.n[1], carry) + t2, carry = bits.Add64(a.n[2], b.n[2], carry) + t3, overflow = bits.Add64(a.n[3], b.n[3], carry) + + // Subtract N with borrow propagation. borrow is set iff t = a+b < N. + // + // u = t - N = a+b - N + var u0, u1, u2, u3, borrow uint64 + u0, borrow = bits.Sub64(t0, orderLimb0, 0) + u1, borrow = bits.Sub64(t1, orderLimb1, borrow) + u2, borrow = bits.Sub64(t2, orderLimb2, borrow) + u3, borrow = bits.Sub64(t3, orderLimb3, borrow) + + // Constant-time select. + // + // Set s = t = a+b only when there was no overflow and t < N (borrow set). + // Otherwise s = u = a+b - N. + cond := (1 - overflow) & borrow + s.n[0] = constantTimeSelect64(cond, t0, u0) + s.n[1] = constantTimeSelect64(cond, t1, u1) + s.n[2] = constantTimeSelect64(cond, t2, u2) + s.n[3] = constantTimeSelect64(cond, t3, u3) return s } @@ -464,315 +384,359 @@ func (s *ModNScalar) Add(val *ModNScalar) *ModNScalar { return s.Add2(s, val) } -// accumulator96 provides a 96-bit accumulator for use in the intermediate -// calculations requiring more than 64-bits. -type accumulator96 struct { - n [3]uint32 -} +// scalar64Reduce512 reduces a 512-bit little-endian limb array modulo the group +// order in constant time and stores the result in r. +func scalar64Reduce512(r *[4]uint64, x *[8]uint64) { + // The overall strategy employed here is: + // 1) Start with the full unreduced 512-bit product of the two scalars. + // 2) Reduce the result modulo the group order via Crandall reduction + // (described below). + // 3) Repeat step 2 noting that each iteration reduces the required number + // of bits by 127 because the two's complement of N has 127 leading zero + // bits. + // 4) Once reduced to less than a max of 2N, perform the final reduction + // with a constant-time conditional subtraction. + // + // Technically the max possible input value here will be (N-1)^2, since the + // only time it is called is with the product of two scalars that are always + // mod N. Nevertheless, it is safer to treat it as the product of two max + // size 256-bit values, (2^256-1)^2 = 2^512 - 2^257 + 1, to ensure + // correctness if that were to ever change. + // + // Per [HAC] section 14.3.4: Reduction method of moduli of special form, + // when the modulus is of the special form m = b^t - c, highly efficient + // reduction can be achieved. While [HAC] only presents the algorithm and + // does not call it out by name or provide the mathematical justification, + // the underlying technique is known as Crandall reduction and is often + // presented as 2^k - c. It is easy to see they are equivalent by setting + // b = 2 and t = k. + // + // The secp256k1 group order is 2^256 - 0x14551231950b75fc4402da1732fc9bebf, + // so it fits this criteria where: + // k = 256 + // c = 432420386565659656852420866394968145599 + // + // Crandall reduction works by taking advantage of the fact that if a prime + // is of the form 2^k - c, then 2^k - c ≡ 0 (mod p), so 2^k ≡ c (mod p). In + // other words, every multiple of 2^k is equivalent to adding c when working + // modulo p. + // + // Since the 512-bit value to reduce is tightly packed into uint64s, the + // upper 4 limbs are all multiples of 2^256. Therefore, reducing modulo the + // group order is equivalent to multiplying those upper limbs by c and + // adding the result to the corresponding lower 4 limbs while propagating + // the carries. + // + // For the specific case of the secp256k1 group order, a max of 4 reductions + // are required because c is 129 bits and so the first round will reduce + // from 512 bits to a max of 385 bits, the second round will reduce to a max + // of 258 bits, and the third round will reduce to within 2N. Then, a + // conditional subtraction of N handles the final reduction. + // + // The reduction is slightly complicated by the fact c needs 3 uint64 limbs + // to represent it, so multiplying the upper limbs by c requires multiplying + // each one by all 3 limbs of c and carefully managing the columns and + // carries. The uppermost limb of c = 1, so multiplication by that limb is + // avoided in the code below. -// Add adds the passed unsigned 64-bit value to the accumulator. -func (a *accumulator96) Add(v uint64) { - low := uint32(v & uint32Mask) - hi := uint32(v >> 32) - a.n[0] += low - hi += constantTimeLess(a.n[0], low) // Carry if overflow in n[0]. - a.n[1] += hi - a.n[2] += constantTimeLess(a.n[1], hi) // Carry if overflow in n[1]. -} + var t0, t1, t2, t3, t4, t5, t6, l0, h0, l1, h1, l2, carry uint64 -// Rsh32 right shifts the accumulator by 32 bits. -func (a *accumulator96) Rsh32() { - a.n[0] = a.n[1] - a.n[1] = a.n[2] - a.n[2] = 0 -} + // ------------------------------------------------------------------------- + // First reduction. + // + // The code in this section is equivalent to the following if uint512 and + // larger numbers and shifts were supported directly: + // + // t0 = uint64(x%(1<<256) + (x>>256)*c) + // t1 = uint64((x%(1<<256) + (x>>256)*c) >> 64) + // t2 = uint64((x%(1<<256) + (x>>256)*c) >> 128) + // t3 = uint64((x%(1<<256) + (x>>256)*c) >> 192) + // t4 = uint64((x%(1<<256) + (x>>256)*c) >> 256) + // t5 = uint64((x%(1<<256) + (x>>256)*c) >> 320) + // t6 = uint64((x%(1<<256) + (x>>256)*c) >> 384) + // + // To visualize the logic, multiply the 3 limbs of c by the upper limbs via + // the standard pencil-and-paper method: + // + // 2^576 2^512 2^448 2^384 2^320 2^256 (place values) + // ----------------------------------- + // x7 x6 x5 x4 + // c2 c1 c0 * + // ----------------------------------- + // x4c2 x4c1 x4c0 + // x5c2 x5c1 x5c0 + + // x6c2 x6c1 x6c0 + + // x7c2 x7c1 x7c0 + + // ----------------------------------- + // x7c2 .... .... .... .... x4c0 + // + // Then add them to the associated lower limbs per the reduction identity: + // + // 2^384 2^320 2^256 2^192 2^128 2^64 2^0 (place values) + // ----------------------------------------- + // x3 x2 x1 x0 + // x4c2 x4c1 x4c0 + + // x5c2 x5c1 x5c0 + + // x6c2 x6c1 x6c0 + + // x7c2 x7c1 x7c0 + + // ----------------------------------------- + // t6 t5 t4 t3 t2 t1 t0 + // + // Where t6 ≤ 1 for the potential carry. + // + // Thus 0 ≤ t < 2^385 after the 1st reduction. + // ------------------------------------------------------------------------- + + // Compute x4*c with result in t0..t2 with carries propagated and the final + // carry in t3. + h0, t0 = bits.Mul64(x[4], orderComplementLimb0) + h1, t1 = bits.Mul64(x[4], orderComplementLimb1) + t1, carry = bits.Add64(t1, h0, 0) + t2, carry = bits.Add64(x[4], h1, carry) + t3 = carry + + // Compute x5*c with result added to t1..t3 with carries propagated and the + // final carry in t4. + // + // Note that since h1 is the upper 64 bits of the product of a uint64 with + // c1 (limb 1 of c) and c1 < 2^63 - 1: + // h1 ≤ floor((2^64-1)(2^63-1) / 2^64) = 2^63 - 2 + // + // Then, because current t3 ≤ 1 and carry ≤ 1, a loose bound for the first + // t3 below is: + // t3 ≤ h1 + 2 = 2^63 < 2^64 + // + // Therefore, it is safe to discard the carry. + h0, l0 = bits.Mul64(x[5], orderComplementLimb0) + h1, l1 = bits.Mul64(x[5], orderComplementLimb1) + t1, carry = bits.Add64(t1, l0, 0) + t2, carry = bits.Add64(t2, h0, carry) + t3, _ = bits.Add64(t3, h1, carry) + t2, carry = bits.Add64(t2, l1, 0) + t3, carry = bits.Add64(t3, x[5], carry) + t4 = carry + + // Compute x6*c with result added to t2..t4 with carries propagated and the + // final carry in t5. + // + // It is safe to discard the carry on the first t4 for the same reason as + // above. + h0, l0 = bits.Mul64(x[6], orderComplementLimb0) + h1, l1 = bits.Mul64(x[6], orderComplementLimb1) + t2, carry = bits.Add64(t2, l0, 0) + t3, carry = bits.Add64(t3, h0, carry) + t4, _ = bits.Add64(t4, h1, carry) + t3, carry = bits.Add64(t3, l1, 0) + t4, carry = bits.Add64(t4, x[6], carry) + t5 = carry + + // Compute x7*c with result added to t3..t5 with carries propagated and the + // final carry in t6. + // + // It is safe to discard the carry on the first t5 for the same reason as + // above. + h0, l0 = bits.Mul64(x[7], orderComplementLimb0) + h1, l1 = bits.Mul64(x[7], orderComplementLimb1) + t3, carry = bits.Add64(t3, l0, 0) + t4, carry = bits.Add64(t4, h0, carry) + t5, _ = bits.Add64(t5, h1, carry) + t4, carry = bits.Add64(t4, l1, 0) + t5, carry = bits.Add64(t5, x[7], carry) + t6 = carry + + // Add result to lower limbs and propagate carries. + // + // It is safe to discard the carry on t6 since the resulting t6 ≤ 1 as + // previously described. + t0, carry = bits.Add64(t0, x[0], 0) + t1, carry = bits.Add64(t1, x[1], carry) + t2, carry = bits.Add64(t2, x[2], carry) + t3, carry = bits.Add64(t3, x[3], carry) + t4, carry = bits.Add64(t4, 0, carry) + t5, carry = bits.Add64(t5, 0, carry) + t6, _ = bits.Add64(t6, 0, carry) + + // ------------------------------------------------------------------------- + // Second reduction. + // + // The value now fits in 385 bits, so reduce it again. Only t4, t5, and t6 + // need to be considered since the higher limb, t7, is ≥ 448 bits and thus + // guaranteed to be 0. Further, as previously noted, t6 ≤ 1. + // + // The code in this section is equivalent to following if uint512 and larger + // numbers and shifts were supported: + // + // t0 = uint64(t%2^256 + (t>>256)*c) + // t1 = uint64((t%2^256 + (t>>256)*c) >> 64) + // t2 = uint64((t%2^256 + (t>>256)*c) >> 128) + // t3 = uint64((t%2^256 + (t>>256)*c) >> 192) + // t4 = uint64((t%2^256 + (t>>256)*c) >> 256) + // + // To visualize the logic, multiply the 3 limbs of c by the new upper limbs + // via the standard pencil-and-paper method: + // + // 2^512 2^448 2^384 2^320 2^256 (place values) + // ----------------------------- + // t6 t5 t4 + // c2 c1 c0 * + // ----------------------------- + // t4c2 t4c1 t4c0 + // t5c2 t5c1 t5c0 + + // t6c2 t6c1 t6c0 + + // ---------------------------- + // t6c2 .... .... .... t4c0 + // + // Then add them to the associated lower limbs per the reduction identity + // noting that the code reuses t for the sums: + // + // 2^256 2^192 2^128 2^64 2^0 (place values) + // ----------------------------- + // t3 t2 t1 t0 + // t4c2 t4c1 t4c0 + + // t5c2 t5c1 t5c0 + + // t6c2 t6c1 t6c0 + + // ----------------------------- + // t4 t3 t2 t1 t0 + // + // With a loose bound for t4 ≤ 3. + // + // Proof: + // + // Known: c2=1, t6 ≤ 1, t3 ≤ 2^64-1, t5 ≤ 2^64-1, and c1 ≤ 2^63-1. + // + // Let s = t3 + t5c2 + t6c1. Then: + // s ≤ 2^64-1 + (2^64-1)*1 + 1*(2^63-1) = 2^65 + 2^63 - 3 + // + // So, the loose bound on the carry in to t4 is carryIn ≤ floor(s/2^64) = 2. + // + // Finally, t4 ≤ t6*c2 + carryIn ≤ 1*1 + 2 ≤ 3. + // + // Thus 0 ≤ t < 2^258 after the 2nd reduction. + // ------------------------------------------------------------------------- -// reduce385 reduces the 385-bit intermediate result in the passed terms modulo -// the group order in constant time and stores the result in s. -func (s *ModNScalar) reduce385(t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12 uint64) { - // At this point, the intermediate result in the passed terms has been - // reduced to fit within 385 bits, so reduce it again using the same method - // described in reduce512. As before, the intermediate result will end up - // being reduced by another 127 bits to 258 bits, thus 9 32-bit terms are - // needed for this iteration. The reduced terms are assigned back to t0 - // through t8. - // - // Note that several of the intermediate calculations require adding 64-bit - // products together which would overflow a uint64, so a 96-bit accumulator - // is used instead until the value is reduced enough to use native uint64s. - - // Terms for 2^(32*0). - var acc accumulator96 - acc.n[0] = uint32(t0) // == acc.Add(t0) because acc is guaranteed to be 0. - acc.Add(t8 * uint64(orderComplementWordZero)) - t0 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*1). - acc.Add(t1) - acc.Add(t8 * uint64(orderComplementWordOne)) - acc.Add(t9 * uint64(orderComplementWordZero)) - t1 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*2). - acc.Add(t2) - acc.Add(t8 * uint64(orderComplementWordTwo)) - acc.Add(t9 * uint64(orderComplementWordOne)) - acc.Add(t10 * uint64(orderComplementWordZero)) - t2 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*3). - acc.Add(t3) - acc.Add(t8 * uint64(orderComplementWordThree)) - acc.Add(t9 * uint64(orderComplementWordTwo)) - acc.Add(t10 * uint64(orderComplementWordOne)) - acc.Add(t11 * uint64(orderComplementWordZero)) - t3 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*4). - acc.Add(t4) - acc.Add(t8) // * uint64(orderComplementWordFour) // * 1 - acc.Add(t9 * uint64(orderComplementWordThree)) - acc.Add(t10 * uint64(orderComplementWordTwo)) - acc.Add(t11 * uint64(orderComplementWordOne)) - acc.Add(t12 * uint64(orderComplementWordZero)) - t4 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*5). - acc.Add(t5) - // acc.Add(t8 * uint64(orderComplementWordFive)) // 0 - acc.Add(t9) // * uint64(orderComplementWordFour) // * 1 - acc.Add(t10 * uint64(orderComplementWordThree)) - acc.Add(t11 * uint64(orderComplementWordTwo)) - acc.Add(t12 * uint64(orderComplementWordOne)) - t5 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*6). - acc.Add(t6) - // acc.Add(t8 * uint64(orderComplementWordSix)) // 0 - // acc.Add(t9 * uint64(orderComplementWordFive)) // 0 - acc.Add(t10) // * uint64(orderComplementWordFour) // * 1 - acc.Add(t11 * uint64(orderComplementWordThree)) - acc.Add(t12 * uint64(orderComplementWordTwo)) - t6 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*7). - acc.Add(t7) - // acc.Add(t8 * uint64(orderComplementWordSeven)) // 0 - // acc.Add(t9 * uint64(orderComplementWordSix)) // 0 - // acc.Add(t10 * uint64(orderComplementWordFive)) // 0 - acc.Add(t11) // * uint64(orderComplementWordFour) // * 1 - acc.Add(t12 * uint64(orderComplementWordThree)) - t7 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*8). - // acc.Add(t9 * uint64(orderComplementWordSeven)) // 0 - // acc.Add(t10 * uint64(orderComplementWordSix)) // 0 - // acc.Add(t11 * uint64(orderComplementWordFive)) // 0 - acc.Add(t12) // * uint64(orderComplementWordFour) // * 1 - t8 = uint64(acc.n[0]) - // acc.Rsh32() // No need since not used after this. Guaranteed to be 0. - - // NOTE: All of the remaining multiplications for this iteration result in 0 - // as they all involve multiplying by combinations of the fifth, sixth, and - // seventh words of the two's complement of N, which are 0, so skip them. - - // At this point, the result is reduced to fit within 258 bits, so reduce it - // again using a slightly modified version of the same method. The maximum - // value in t8 is 2 at this point and therefore multiplying it by each word - // of the two's complement of N and adding it to a 32-bit term will result - // in a maximum requirement of 33 bits, so it is safe to use native uint64s - // here for the intermediate term carry propagation. - // - // Also, since the maximum value in t8 is 2, this ends up reducing by - // another 2 bits to 256 bits. - c := t0 + t8*uint64(orderComplementWordZero) - s.n[0] = uint32(c & uint32Mask) - c = (c >> 32) + t1 + t8*uint64(orderComplementWordOne) - s.n[1] = uint32(c & uint32Mask) - c = (c >> 32) + t2 + t8*uint64(orderComplementWordTwo) - s.n[2] = uint32(c & uint32Mask) - c = (c >> 32) + t3 + t8*uint64(orderComplementWordThree) - s.n[3] = uint32(c & uint32Mask) - c = (c >> 32) + t4 + t8 // * uint64(orderComplementWordFour) == * 1 - s.n[4] = uint32(c & uint32Mask) - c = (c >> 32) + t5 // + t8*uint64(orderComplementWordFive) == 0 - s.n[5] = uint32(c & uint32Mask) - c = (c >> 32) + t6 // + t8*uint64(orderComplementWordSix) == 0 - s.n[6] = uint32(c & uint32Mask) - c = (c >> 32) + t7 // + t8*uint64(orderComplementWordSeven) == 0 - s.n[7] = uint32(c & uint32Mask) - - // The result is now 256 bits, but it might still be >= N, so use the - // existing normal reduce method for 256-bit values. - s.reduce256(uint32(c>>32) + s.overflows()) -} + // Compute t4*c with result added to t0..t2 with carries propagated and the + // final carry in t4. + // + // It is safe to discard the carry on the final t4 given the carries ≤ 1. + h0, l0 = bits.Mul64(t4, orderComplementLimb0) + h1, l1 = bits.Mul64(t4, orderComplementLimb1) + l2 = t4 // h2, l2 = bits.Mul64(t4, orderComplementLimb2) => h2=0, l2=t4 + t0, carry = bits.Add64(t0, l0, 0) + t1, carry = bits.Add64(t1, h0, carry) + t2, carry = bits.Add64(t2, h1, carry) + t3, carry = bits.Add64(t3, 0, carry) + t4 = carry + t1, carry = bits.Add64(t1, l1, 0) + t2, carry = bits.Add64(t2, l2, carry) + t3, carry = bits.Add64(t3, 0, carry) + t4, _ = bits.Add64(t4, 0, carry) + + // Compute t5*c with result added to t1..t3 with carries propagated. + // + // It is safe to discard the carry on the final t4 given the carries ≤ 1. + h0, l0 = bits.Mul64(t5, orderComplementLimb0) + h1, l1 = bits.Mul64(t5, orderComplementLimb1) + // h2, l2 = bits.Mul64(t5, orderComplementLimb2) => h2=0, l2=t5 + t1, carry = bits.Add64(t1, l0, 0) + t2, carry = bits.Add64(t2, h0, carry) + t3, carry = bits.Add64(t3, h1, carry) + t4, _ = bits.Add64(t4, 0, carry) + t2, carry = bits.Add64(t2, l1, 0) + t3, carry = bits.Add64(t3, t5, carry) + t4, _ = bits.Add64(t4, 0, carry) + + // Compute t6*c with result added to t2..t4 with carries propagated. + // + // Note that since t6 ≤ 1, the product can't overflow a uint64, so it is + // safe to use normal 64-bit multiplication. + // + // It is safe to discard the carry on t4 since the resulting t4 ≤ 3 as + // previously described. + t2, carry = bits.Add64(t2, t6*orderComplementLimb0, 0) + t3, carry = bits.Add64(t3, t6*orderComplementLimb1, carry) + t4, _ = bits.Add64(t4, t6, carry) + + // ------------------------------------------------------------------------- + // Third reduction. + // + // The value now fits in 258 bits, so reduce it again. Only t4 needs to be + // considered since the higher limbs are ≥ 320 bits and thus guaranteed to + // be 0. + // + // The code in this section is equivalent to following if uint512 and larger + // numbers and shifts were supported: + // + // t0 = uint64(t%2^256 + (t>>256)*c) + // t1 = uint64((t%2^256 + (t>>256)*c) >> 64) + // t2 = uint64((t%2^256 + (t>>256)*c) >> 128) + // t3 = uint64((t%2^256 + (t>>256)*c) >> 192) + // t4 = uint64((t%2^256 + (t>>256)*c) >> 256) + // + // To visualize the logic, multiply the 3 limbs of c by the new upper limb + // via the standard pencil-and-paper method: + // + // 2^512 2^448 2^384 2^320 2^256 (place values) + // ----------------------------- + // t4 + // c2 c1 c0 * + // ----------------------------- + // t4c2 t4c1 t4c0 + // + // Then add them to the associated lower limbs noting that the code reuses t + // for the sums: + // + // 2^256 2^192 2^128 2^64 2^0 (place values) + // ----------------------------- + // t3 t2 t1 t0 + // t4c2 t4c1 t4c0 + + // ----------------------------- + // t4 t3 t2 t1 t0 + // + // Where t4 ≤ 1 for the carry. + // ------------------------------------------------------------------------- -// reduce512 reduces the 512-bit intermediate result in the passed terms modulo -// the group order down to 385 bits in constant time and stores the result in s. -func (s *ModNScalar) reduce512(t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15 uint64) { - // At this point, the intermediate result in the passed terms is grouped - // into the respective bases. + // Compute t4*c with result added to t0..t2 with carries propagated and the + // final carry in t4. // - // Per [HAC] section 14.3.4: Reduction method of moduli of special form, - // when the modulus is of the special form m = b^t - c, where log_2(c) < t, - // highly efficient reduction can be achieved per the provided algorithm. - // - // The secp256k1 group order fits this criteria since it is: - // 2^256 - 432420386565659656852420866394968145599 - // - // Technically the max possible value here is (N-1)^2 since the two scalars - // being multiplied are always mod N. Nevertheless, it is safer to consider - // it to be (2^256-1)^2 = 2^512 - 2^257 + 1 since it is the product of two - // 256-bit values. - // - // The algorithm is to reduce the result modulo the prime by subtracting - // multiples of the group order N. However, in order simplify carry - // propagation, this adds with the two's complement of N to achieve the same - // result. - // - // Since the two's complement of N has 127 leading zero bits, this will end - // up reducing the intermediate result from 512 bits to 385 bits, resulting - // in 13 32-bit terms. The reduced terms are assigned back to t0 through - // t12. - // - // Note that several of the intermediate calculations require adding 64-bit - // products together which would overflow a uint64, so a 96-bit accumulator - // is used instead. - - // Terms for 2^(32*0). - var acc accumulator96 - acc.n[0] = uint32(t0) // == acc.Add(t0) because acc is guaranteed to be 0. - acc.Add(t8 * uint64(orderComplementWordZero)) - t0 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*1). - acc.Add(t1) - acc.Add(t8 * uint64(orderComplementWordOne)) - acc.Add(t9 * uint64(orderComplementWordZero)) - t1 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*2). - acc.Add(t2) - acc.Add(t8 * uint64(orderComplementWordTwo)) - acc.Add(t9 * uint64(orderComplementWordOne)) - acc.Add(t10 * uint64(orderComplementWordZero)) - t2 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*3). - acc.Add(t3) - acc.Add(t8 * uint64(orderComplementWordThree)) - acc.Add(t9 * uint64(orderComplementWordTwo)) - acc.Add(t10 * uint64(orderComplementWordOne)) - acc.Add(t11 * uint64(orderComplementWordZero)) - t3 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*4). - acc.Add(t4) - acc.Add(t8) // * uint64(orderComplementWordFour) // * 1 - acc.Add(t9 * uint64(orderComplementWordThree)) - acc.Add(t10 * uint64(orderComplementWordTwo)) - acc.Add(t11 * uint64(orderComplementWordOne)) - acc.Add(t12 * uint64(orderComplementWordZero)) - t4 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*5). - acc.Add(t5) - // acc.Add(t8 * uint64(orderComplementWordFive)) // 0 - acc.Add(t9) // * uint64(orderComplementWordFour) // * 1 - acc.Add(t10 * uint64(orderComplementWordThree)) - acc.Add(t11 * uint64(orderComplementWordTwo)) - acc.Add(t12 * uint64(orderComplementWordOne)) - acc.Add(t13 * uint64(orderComplementWordZero)) - t5 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*6). - acc.Add(t6) - // acc.Add(t8 * uint64(orderComplementWordSix)) // 0 - // acc.Add(t9 * uint64(orderComplementWordFive)) // 0 - acc.Add(t10) // * uint64(orderComplementWordFour)) // * 1 - acc.Add(t11 * uint64(orderComplementWordThree)) - acc.Add(t12 * uint64(orderComplementWordTwo)) - acc.Add(t13 * uint64(orderComplementWordOne)) - acc.Add(t14 * uint64(orderComplementWordZero)) - t6 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*7). - acc.Add(t7) - // acc.Add(t8 * uint64(orderComplementWordSeven)) // 0 - // acc.Add(t9 * uint64(orderComplementWordSix)) // 0 - // acc.Add(t10 * uint64(orderComplementWordFive)) // 0 - acc.Add(t11) // * uint64(orderComplementWordFour) // * 1 - acc.Add(t12 * uint64(orderComplementWordThree)) - acc.Add(t13 * uint64(orderComplementWordTwo)) - acc.Add(t14 * uint64(orderComplementWordOne)) - acc.Add(t15 * uint64(orderComplementWordZero)) - t7 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*8). - // acc.Add(t9 * uint64(orderComplementWordSeven)) // 0 - // acc.Add(t10 * uint64(orderComplementWordSix)) // 0 - // acc.Add(t11 * uint64(orderComplementWordFive)) // 0 - acc.Add(t12) // * uint64(orderComplementWordFour) // * 1 - acc.Add(t13 * uint64(orderComplementWordThree)) - acc.Add(t14 * uint64(orderComplementWordTwo)) - acc.Add(t15 * uint64(orderComplementWordOne)) - t8 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*9). - // acc.Add(t10 * uint64(orderComplementWordSeven)) // 0 - // acc.Add(t11 * uint64(orderComplementWordSix)) // 0 - // acc.Add(t12 * uint64(orderComplementWordFive)) // 0 - acc.Add(t13) // * uint64(orderComplementWordFour) // * 1 - acc.Add(t14 * uint64(orderComplementWordThree)) - acc.Add(t15 * uint64(orderComplementWordTwo)) - t9 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*10). - // acc.Add(t11 * uint64(orderComplementWordSeven)) // 0 - // acc.Add(t12 * uint64(orderComplementWordSix)) // 0 - // acc.Add(t13 * uint64(orderComplementWordFive)) // 0 - acc.Add(t14) // * uint64(orderComplementWordFour) // * 1 - acc.Add(t15 * uint64(orderComplementWordThree)) - t10 = uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*11). - // acc.Add(t12 * uint64(orderComplementWordSeven)) // 0 - // acc.Add(t13 * uint64(orderComplementWordSix)) // 0 - // acc.Add(t14 * uint64(orderComplementWordFive)) // 0 - acc.Add(t15) // * uint64(orderComplementWordFour) // * 1 - t11 = uint64(acc.n[0]) - acc.Rsh32() - - // NOTE: All of the remaining multiplications for this iteration result in 0 - // as they all involve multiplying by combinations of the fifth, sixth, and - // seventh words of the two's complement of N, which are 0, so skip them. - - // Terms for 2^(32*12). - t12 = uint64(acc.n[0]) - // acc.Rsh32() // No need since not used after this. Guaranteed to be 0. - - // At this point, the result is reduced to fit within 385 bits, so reduce it - // again using the same method accordingly. - s.reduce385(t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12) + // Note that since current t4 ≤ 3, these bounds for the product hold: + // t4*c0 < 2^64 + // t4*c1 < 2^64 + // + // So, uint64 overflow is impossible and therefore it is safe to use normal + // 64-bit multiplication. + t0, carry = bits.Add64(t0, t4*orderComplementLimb0, 0) + t1, carry = bits.Add64(t1, t4*orderComplementLimb1, carry) + t2, carry = bits.Add64(t2, t4, carry) + t3, carry = bits.Add64(t3, 0, carry) + t4 = carry + + // ------------------------------------------------------------------------- + // Final reduction. + // + // The value is now in the range 0 ≤ t < 2N, so one 5-limb conditional + // subtract of N guarantees it is fully reduced. + // ------------------------------------------------------------------------- + + // Subtract N with borrow propagation. borrow is set iff t < N. + // + // In other words, the value overflowed (≥ N) when t - N does NOT borrow. + // + // s = t - N + var s0, s1, s2, s3, borrow uint64 + s0, borrow = bits.Sub64(t0, orderLimb0, 0) + s1, borrow = bits.Sub64(t1, orderLimb1, borrow) + s2, borrow = bits.Sub64(t2, orderLimb2, borrow) + s3, borrow = bits.Sub64(t3, orderLimb3, borrow) + _, borrow = bits.Sub64(t4, 0, borrow) + + // Constant-time select. + // + // Set r = t only when t < N (borrow set). + // Otherwise r = s = t - N. + r[0] = constantTimeSelect64(borrow, t0, s0) + r[1] = constantTimeSelect64(borrow, t1, s1) + r[2] = constantTimeSelect64(borrow, t2, s2) + r[3] = constantTimeSelect64(borrow, t3, s3) } // Mul2 multiplies the passed two scalars together modulo the group order in @@ -780,159 +744,10 @@ func (s *ModNScalar) reduce512(t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, // // The scalar is returned to support chaining. This enables syntax like: // s3.Mul2(s, s2).AddInt(1) so that s3 = (s * s2) + 1. -func (s *ModNScalar) Mul2(val, val2 *ModNScalar) *ModNScalar { - // This could be done with for loops and an array to store the intermediate - // terms, but this unrolled version is significantly faster. - - // The overall strategy employed here is: - // 1) Calculate the 512-bit product of the two scalars using the standard - // pencil-and-paper method. - // 2) Reduce the result modulo the prime by effectively subtracting - // multiples of the group order N (actually performed by adding multiples - // of the two's complement of N to avoid implementing subtraction). - // 3) Repeat step 2 noting that each iteration reduces the required number - // of bits by 127 because the two's complement of N has 127 leading zero - // bits. - // 4) Once reduced to 256 bits, call the existing reduce method to perform - // a final reduction as needed. - // - // Note that several of the intermediate calculations require adding 64-bit - // products together which would overflow a uint64, so a 96-bit accumulator - // is used instead. - - // Terms for 2^(32*0). - var acc accumulator96 - acc.Add(uint64(val.n[0]) * uint64(val2.n[0])) - t0 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*1). - acc.Add(uint64(val.n[0]) * uint64(val2.n[1])) - acc.Add(uint64(val.n[1]) * uint64(val2.n[0])) - t1 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*2). - acc.Add(uint64(val.n[0]) * uint64(val2.n[2])) - acc.Add(uint64(val.n[1]) * uint64(val2.n[1])) - acc.Add(uint64(val.n[2]) * uint64(val2.n[0])) - t2 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*3). - acc.Add(uint64(val.n[0]) * uint64(val2.n[3])) - acc.Add(uint64(val.n[1]) * uint64(val2.n[2])) - acc.Add(uint64(val.n[2]) * uint64(val2.n[1])) - acc.Add(uint64(val.n[3]) * uint64(val2.n[0])) - t3 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*4). - acc.Add(uint64(val.n[0]) * uint64(val2.n[4])) - acc.Add(uint64(val.n[1]) * uint64(val2.n[3])) - acc.Add(uint64(val.n[2]) * uint64(val2.n[2])) - acc.Add(uint64(val.n[3]) * uint64(val2.n[1])) - acc.Add(uint64(val.n[4]) * uint64(val2.n[0])) - t4 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*5). - acc.Add(uint64(val.n[0]) * uint64(val2.n[5])) - acc.Add(uint64(val.n[1]) * uint64(val2.n[4])) - acc.Add(uint64(val.n[2]) * uint64(val2.n[3])) - acc.Add(uint64(val.n[3]) * uint64(val2.n[2])) - acc.Add(uint64(val.n[4]) * uint64(val2.n[1])) - acc.Add(uint64(val.n[5]) * uint64(val2.n[0])) - t5 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*6). - acc.Add(uint64(val.n[0]) * uint64(val2.n[6])) - acc.Add(uint64(val.n[1]) * uint64(val2.n[5])) - acc.Add(uint64(val.n[2]) * uint64(val2.n[4])) - acc.Add(uint64(val.n[3]) * uint64(val2.n[3])) - acc.Add(uint64(val.n[4]) * uint64(val2.n[2])) - acc.Add(uint64(val.n[5]) * uint64(val2.n[1])) - acc.Add(uint64(val.n[6]) * uint64(val2.n[0])) - t6 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*7). - acc.Add(uint64(val.n[0]) * uint64(val2.n[7])) - acc.Add(uint64(val.n[1]) * uint64(val2.n[6])) - acc.Add(uint64(val.n[2]) * uint64(val2.n[5])) - acc.Add(uint64(val.n[3]) * uint64(val2.n[4])) - acc.Add(uint64(val.n[4]) * uint64(val2.n[3])) - acc.Add(uint64(val.n[5]) * uint64(val2.n[2])) - acc.Add(uint64(val.n[6]) * uint64(val2.n[1])) - acc.Add(uint64(val.n[7]) * uint64(val2.n[0])) - t7 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*8). - acc.Add(uint64(val.n[1]) * uint64(val2.n[7])) - acc.Add(uint64(val.n[2]) * uint64(val2.n[6])) - acc.Add(uint64(val.n[3]) * uint64(val2.n[5])) - acc.Add(uint64(val.n[4]) * uint64(val2.n[4])) - acc.Add(uint64(val.n[5]) * uint64(val2.n[3])) - acc.Add(uint64(val.n[6]) * uint64(val2.n[2])) - acc.Add(uint64(val.n[7]) * uint64(val2.n[1])) - t8 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*9). - acc.Add(uint64(val.n[2]) * uint64(val2.n[7])) - acc.Add(uint64(val.n[3]) * uint64(val2.n[6])) - acc.Add(uint64(val.n[4]) * uint64(val2.n[5])) - acc.Add(uint64(val.n[5]) * uint64(val2.n[4])) - acc.Add(uint64(val.n[6]) * uint64(val2.n[3])) - acc.Add(uint64(val.n[7]) * uint64(val2.n[2])) - t9 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*10). - acc.Add(uint64(val.n[3]) * uint64(val2.n[7])) - acc.Add(uint64(val.n[4]) * uint64(val2.n[6])) - acc.Add(uint64(val.n[5]) * uint64(val2.n[5])) - acc.Add(uint64(val.n[6]) * uint64(val2.n[4])) - acc.Add(uint64(val.n[7]) * uint64(val2.n[3])) - t10 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*11). - acc.Add(uint64(val.n[4]) * uint64(val2.n[7])) - acc.Add(uint64(val.n[5]) * uint64(val2.n[6])) - acc.Add(uint64(val.n[6]) * uint64(val2.n[5])) - acc.Add(uint64(val.n[7]) * uint64(val2.n[4])) - t11 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*12). - acc.Add(uint64(val.n[5]) * uint64(val2.n[7])) - acc.Add(uint64(val.n[6]) * uint64(val2.n[6])) - acc.Add(uint64(val.n[7]) * uint64(val2.n[5])) - t12 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*13). - acc.Add(uint64(val.n[6]) * uint64(val2.n[7])) - acc.Add(uint64(val.n[7]) * uint64(val2.n[6])) - t13 := uint64(acc.n[0]) - acc.Rsh32() - - // Terms for 2^(32*14). - acc.Add(uint64(val.n[7]) * uint64(val2.n[7])) - t14 := uint64(acc.n[0]) - acc.Rsh32() - - // What's left is for 2^(32*15). - t15 := uint64(acc.n[0]) - // acc.Rsh32() // No need since not used after this. Guaranteed to be 0. - - // At this point, all of the terms are grouped into their respective base - // and occupy up to 512 bits. Reduce the result accordingly. - s.reduce512(t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, - t15) +func (s *ModNScalar) Mul2(a, b *ModNScalar) *ModNScalar { + var product [8]uint64 + field64Mul512(&product, &a.n, &b.n) + scalar64Reduce512(&s.n, &product) return s } @@ -951,15 +766,10 @@ func (s *ModNScalar) Mul(val *ModNScalar) *ModNScalar { // The scalar is returned to support chaining. This enables syntax like: // s3.SquareVal(s).Mul(s) so that s3 = s^2 * s = s^3. func (s *ModNScalar) SquareVal(val *ModNScalar) *ModNScalar { - // This could technically be optimized slightly to take advantage of the - // fact that many of the intermediate calculations in squaring are just - // doubling, however, benchmarking has shown that due to the need to use a - // 96-bit accumulator, any savings are essentially offset by that and - // consequently there is no real difference in performance over just - // multiplying the value by itself to justify the extra code for now. This - // can be revisited in the future if it becomes a bottleneck in practice. - - return s.Mul2(val, val) + var product [8]uint64 + field64Square512(&product, &val.n) + scalar64Reduce512(&s.n, &product) + return s } // Square squares the scalar modulo the group order in constant time. The @@ -982,32 +792,38 @@ func (s *ModNScalar) NegateVal(val *ModNScalar) *ModNScalar { // minus the value. This implies that the result will always be in the // desired range with the sole exception of 0 because N - 0 = N itself. // - // Therefore, in order to avoid the need to reduce the result for every - // other case in order to achieve constant time, this creates a mask that is - // all 0s in the case of the scalar being negated is 0 and all 1s otherwise - // and bitwise ands that mask with each word. - // - // Finally, to simplify the carry propagation, this adds the two's - // complement of the scalar to N in order to achieve the same result. - bits := val.n[0] | val.n[1] | val.n[2] | val.n[3] | val.n[4] | val.n[5] | - val.n[6] | val.n[7] - mask := uint64(uint32Mask * constantTimeNotEq(bits, 0)) - c := uint64(orderWordZero) + (uint64(^val.n[0]) + 1) - s.n[0] = uint32(c & mask) - c = (c >> 32) + uint64(orderWordOne) + uint64(^val.n[1]) - s.n[1] = uint32(c & mask) - c = (c >> 32) + uint64(orderWordTwo) + uint64(^val.n[2]) - s.n[2] = uint32(c & mask) - c = (c >> 32) + uint64(orderWordThree) + uint64(^val.n[3]) - s.n[3] = uint32(c & mask) - c = (c >> 32) + uint64(orderWordFour) + uint64(^val.n[4]) - s.n[4] = uint32(c & mask) - c = (c >> 32) + uint64(orderWordFive) + uint64(^val.n[5]) - s.n[5] = uint32(c & mask) - c = (c >> 32) + uint64(orderWordSix) + uint64(^val.n[6]) - s.n[6] = uint32(c & mask) - c = (c >> 32) + uint64(orderWordSeven) + uint64(^val.n[7]) - s.n[7] = uint32(c & mask) + // The following handles that case in constant time by creating a mask that + // is all 0s in the case the scalar being negated is 0 and all 1s otherwise + // and then bitwise ands that mask with each limb. + // + // This approach was chosen over subtracting from 0 and then conditionally + // adding N because it's over twice as fast on 32-bit hardware while only + // being about 3-4% slower on 64-bit hardware. + // + // Determine mask first to allow aliasing. + mask := -uint64(constantTimeNotEq64(val.n[0]|val.n[1]|val.n[2]|val.n[3], 0)) + + // Unconditionally subtract the scalar from the group order. + // + // To simplify the carry propagation, this adds the two's complement of the + // scalar to N in order to achieve the same result. + // + // s = N - val + var c uint64 + s.n[0], c = bits.Add64(orderLimb0, ^val.n[0], 1) + s.n[1], c = bits.Add64(orderLimb1, ^val.n[1], c) + s.n[2], c = bits.Add64(orderLimb2, ^val.n[2], c) + s.n[3], _ = bits.Add64(orderLimb3, ^val.n[3], c) + + // Either keep the result when val != 0 or clear it when val == 0. The + // result is either: + // + // val == 0: s = 0 + // val != 0: s = N - val + s.n[0] &= mask + s.n[1] &= mask + s.n[2] &= mask + s.n[3] &= mask return s } @@ -1079,26 +895,42 @@ func (s *ModNScalar) InverseNonConst() *ModNScalar { // IsOverHalfOrder returns whether or not the scalar exceeds the group order // divided by 2 in constant time. func (s *ModNScalar) IsOverHalfOrder() bool { - // The intuition here is that the scalar is greater than half of the group - // order if one of the higher individual words is greater than the - // corresponding word of the half group order and all higher words in the - // scalar are equal to their corresponding word of the half group order. - // - // Note that the words 4, 5, and 6 are all the max uint32 value, so there is - // no need to test if those individual words of the scalar exceeds them, - // hence, only equality is checked for them. - result := constantTimeGreater(s.n[7], halfOrderWordSeven) - highWordsEqual := constantTimeEq(s.n[7], halfOrderWordSeven) - highWordsEqual &= constantTimeEq(s.n[6], halfOrderWordSix) - highWordsEqual &= constantTimeEq(s.n[5], halfOrderWordFive) - highWordsEqual &= constantTimeEq(s.n[4], halfOrderWordFour) - result |= highWordsEqual & constantTimeGreater(s.n[3], halfOrderWordThree) - highWordsEqual &= constantTimeEq(s.n[3], halfOrderWordThree) - result |= highWordsEqual & constantTimeGreater(s.n[2], halfOrderWordTwo) - highWordsEqual &= constantTimeEq(s.n[2], halfOrderWordTwo) - result |= highWordsEqual & constantTimeGreater(s.n[1], halfOrderWordOne) - highWordsEqual &= constantTimeEq(s.n[1], halfOrderWordOne) - result |= highWordsEqual & constantTimeGreater(s.n[0], halfOrderWordZero) - - return result != 0 + // These fields provide convenient access to each of the limbs of the + // secp256k1 curve group order N / 2 to improve code readability and avoid + // the need to recalculate them. + // + // The half order of the secp256k1 curve group is: + // 0x7fffffffffffffffffffffffffffffff5d576e7357a4501ddfe92f46681b20a0 + // + // Converting that to field representation (base 2^64) is: + // + // n[0] = 0xdfe92f46681b20a0 + // n[1] = 0x5d576e7357a4501d + // n[2] = 0xffffffffffffffff + // n[3] = 0x7fffffffffffffff + // + // This can be verified with the following test code: + // halfOrder := new(big.Int).Div(curveParams.N, big.NewInt(2)) + // var s ModNScalar + // s.SetByteSlice(halfOrder.Bytes()) + // t.Logf("%x", s.n) + const ( + halfOrderLimb0 uint64 = 0xdfe92f46681b20a0 + halfOrderLimb1 uint64 = 0x5d576e7357a4501d + halfOrderLimb2 uint64 = 0xffffffffffffffff + halfOrderLimb3 uint64 = 0x7fffffffffffffff + ) + + // The goal is to return true when the scalar is greater than half of the + // group order. That is, return true when s > N/2, which is trivially + // rearranged to s - (N/2 + 1) ≥ 0. + // + // In other words, the condition is met iff subtracting (N/2 + 1) from s is + // non-negative (aka there was no borrow). + var borrow uint64 + _, borrow = bits.Sub64(s.n[0], halfOrderLimb0+1, borrow) + _, borrow = bits.Sub64(s.n[1], halfOrderLimb1, borrow) + _, borrow = bits.Sub64(s.n[2], halfOrderLimb2, borrow) + _, borrow = bits.Sub64(s.n[3], halfOrderLimb3, borrow) + return borrow == 0 } diff --git a/dcrec/secp256k1/modnscalar_test.go b/dcrec/secp256k1/modnscalar_test.go index 448378e0c..94e329f23 100644 --- a/dcrec/secp256k1/modnscalar_test.go +++ b/dcrec/secp256k1/modnscalar_test.go @@ -1,4 +1,4 @@ -// Copyright (c) 2020-2023 The Decred developers +// Copyright (c) 2020-2026 The Decred developers // Use of this source code is governed by an ISC // license that can be found in the LICENSE file. @@ -116,23 +116,27 @@ func TestModNScalarSetInt(t *testing.T) { tests := []struct { name string // test description in uint32 // test value - expected [8]uint32 // expected raw ints + expected [4]uint64 // expected limbs }{{ + name: "zero", + in: 0, + expected: [4]uint64{0, 0, 0, 0}, + }, { name: "five", in: 5, - expected: [8]uint32{5, 0, 0, 0, 0, 0, 0, 0}, + expected: [4]uint64{5, 0, 0, 0}, }, { - name: "group order word zero", - in: orderWordZero, - expected: [8]uint32{orderWordZero, 0, 0, 0, 0, 0, 0, 0}, + name: "max uint8 (2^8 - 1)", + in: 255, + expected: [4]uint64{255, 0, 0, 0}, }, { - name: "group order word zero + 1", - in: orderWordZero + 1, - expected: [8]uint32{orderWordZero + 1, 0, 0, 0, 0, 0, 0, 0}, + name: "max uint16 (2^16 - 1)", + in: 65535, + expected: [4]uint64{65535, 0, 0, 0}, }, { - name: "2^32 - 1", + name: "max uint32 (2^32 - 1)", in: 4294967295, - expected: [8]uint32{4294967295, 0, 0, 0, 0, 0, 0, 0}, + expected: [4]uint64{4294967295, 0, 0, 0}, }} for _, test := range tests { @@ -152,88 +156,136 @@ func TestModNScalarSetBytes(t *testing.T) { tests := []struct { name string // test description in string // hex encoded test value - expected [8]uint32 // expected raw ints + expected [4]uint64 // expected raw ints overflow bool // expected overflow result }{{ name: "zero", in: "00", - expected: [8]uint32{0, 0, 0, 0, 0, 0, 0, 0}, + expected: [4]uint64{0, 0, 0, 0}, + overflow: false, + }, { + name: "one (only limb zero set)", + in: "01", + expected: [4]uint64{1, 0, 0, 0}, + overflow: false, + }, { + name: "2^64 - 1 (limb zero all ones)", + in: "ffffffffffffffff", + expected: [4]uint64{0xffffffffffffffff, 0, 0, 0}, + overflow: false, + }, { + name: "2^64 (limb one low bit)", + in: "010000000000000000", + expected: [4]uint64{0, 1, 0, 0}, + overflow: false, + }, { + name: "2^128 - 1 (limbs zero and one all ones)", + in: "ffffffffffffffffffffffffffffffff", + expected: [4]uint64{0xffffffffffffffff, 0xffffffffffffffff, 0, 0}, + overflow: false, + }, { + name: "2^128 (limb two low bit)", + in: "0100000000000000000000000000000000", + expected: [4]uint64{0, 0, 1, 0}, + overflow: false, + }, { + name: "2^192 - 1 (limbs zero, one, and two all ones)", + in: "ffffffffffffffffffffffffffffffffffffffffffffffff", + expected: [4]uint64{ + 0xffffffffffffffff, 0xffffffffffffffff, 0xffffffffffffffff, 0, + }, + overflow: false, + }, { + name: "2^192 (limb three low bit)", + in: "01000000000000000000000000000000000000000000000000", + expected: [4]uint64{0, 0, 0, 1}, + overflow: false, + }, { + name: "2^255 (limb three high bit)", + in: "8000000000000000000000000000000000000000000000000000000000000000", + expected: [4]uint64{0, 0, 0, 0x8000000000000000}, + overflow: false, + }, { + name: "distinct value in every limb (verifies limb ordering)", + in: "0123456789abcdeffedcba98765432100f1e2d3c4b5a69788796a5b4c3d2e1f0", + expected: [4]uint64{ + 0x08796a5b4c3d2e1f0, 0x0f1e2d3c4b5a6978, 0xfedcba9876543210, + 0x0123456789abcdef, + }, overflow: false, }, { name: "group order (aka 0)", in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141", - expected: [8]uint32{0, 0, 0, 0, 0, 0, 0, 0}, + expected: [4]uint64{0, 0, 0, 0}, overflow: true, }, { - name: "group order - 1", + name: "group order - 1 (limb zero one below prime, upper limbs maxed)", in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", - expected: [8]uint32{ - 0xd0364140, 0xbfd25e8c, 0xaf48a03b, 0xbaaedce6, - 0xfffffffe, 0xffffffff, 0xffffffff, 0xffffffff, + expected: [4]uint64{ + 0xbfd25e8cd0364140, 0xbaaedce6af48a03b, 0xfffffffffffffffe, + 0xffffffffffffffff, }, overflow: false, }, { - name: "group order + 1 (aka 1, overflow in word zero)", + name: "group order + 1 (aka 1, overflow in limb zero)", in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364142", - expected: [8]uint32{1, 0, 0, 0, 0, 0, 0, 0}, + expected: [4]uint64{1, 0, 0, 0}, overflow: true, }, { - name: "group order word zero", - in: "d0364141", - expected: [8]uint32{0xd0364141, 0, 0, 0, 0, 0, 0, 0}, - overflow: false, - }, { - name: "group order word zero and one", + name: "group order limb zero", in: "bfd25e8cd0364141", - expected: [8]uint32{0xd0364141, 0xbfd25e8c, 0, 0, 0, 0, 0, 0}, + expected: [4]uint64{0xbfd25e8cd0364141, 0, 0, 0}, overflow: false, }, { - name: "group order words zero, one, and two", - in: "af48a03bbfd25e8cd0364141", - expected: [8]uint32{0xd0364141, 0xbfd25e8c, 0xaf48a03b, 0, 0, 0, 0, 0}, + name: "group order limbs zero and one", + in: "baaedce6af48a03bbfd25e8cd0364141", + expected: [4]uint64{0xbfd25e8cd0364141, 0xbaaedce6af48a03b, 0, 0}, overflow: false, }, { - name: "overflow in word one", - in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8dd0364141", - expected: [8]uint32{0, 1, 0, 0, 0, 0, 0, 0}, - overflow: true, + name: "group order limbs zero, one, and two", + in: "fffffffffffffffebaaedce6af48a03bbfd25e8cd0364141", + expected: [4]uint64{ + 0xbfd25e8cd0364141, 0xbaaedce6af48a03b, 0xfffffffffffffffe, 0, + }, + overflow: false, }, { - name: "overflow in word two", + name: "group order + 2^64 (overflow in limb one)", in: "fffffffffffffffffffffffffffffffebaaedce6af48a03cbfd25e8cd0364141", - expected: [8]uint32{0, 0, 1, 0, 0, 0, 0, 0}, + expected: [4]uint64{0, 1, 0, 0}, overflow: true, }, { - name: "overflow in word three", - in: "fffffffffffffffffffffffffffffffebaaedce7af48a03bbfd25e8cd0364141", - expected: [8]uint32{0, 0, 0, 1, 0, 0, 0, 0}, + name: "group order + 2^128 (overflow in limb two)", + in: "ffffffffffffffffffffffffffffffffbaaedce6af48a03bbfd25e8cd0364141", + expected: [4]uint64{0, 0, 1, 0}, overflow: true, }, { - name: "overflow in word four", - in: "ffffffffffffffffffffffffffffffffbaaedce6af48a03bbfd25e8cd0364141", - expected: [8]uint32{0, 0, 0, 0, 1, 0, 0, 0}, + name: "2^256 - 1 (max input)", + in: "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", + expected: [4]uint64{0x402da1732fc9bebe, 0x4551231950b75fc4, 1, 0}, + overflow: true, }, { name: "(group order - 1) * 2 NOT mod N, truncated >32 bytes", in: "01fffffffffffffffffffffffffffffffd755db9cd5e9140777fa4bd19a06c8284", - expected: [8]uint32{ - 0x19a06c82, 0x777fa4bd, 0xcd5e9140, 0xfd755db9, - 0xffffffff, 0xffffffff, 0xffffffff, 0x01ffffff, + expected: [4]uint64{ + 0x777fa4bd19a06c82, 0xfd755db9cd5e9140, 0xffffffffffffffff, + 0x01ffffffffffffff, }, overflow: false, }, { - name: "alternating bits", + name: "alternating bits (spans all limbs)", in: "a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5", - expected: [8]uint32{ - 0xa5a5a5a5, 0xa5a5a5a5, 0xa5a5a5a5, 0xa5a5a5a5, - 0xa5a5a5a5, 0xa5a5a5a5, 0xa5a5a5a5, 0xa5a5a5a5, + expected: [4]uint64{ + 0xa5a5a5a5a5a5a5a5, 0xa5a5a5a5a5a5a5a5, 0xa5a5a5a5a5a5a5a5, + 0xa5a5a5a5a5a5a5a5, }, overflow: false, }, { - name: "alternating bits 2", + name: "alternating bits 2 (spans all limbs)", in: "5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a", - expected: [8]uint32{ - 0x5a5a5a5a, 0x5a5a5a5a, 0x5a5a5a5a, 0x5a5a5a5a, - 0x5a5a5a5a, 0x5a5a5a5a, 0x5a5a5a5a, 0x5a5a5a5a, + expected: [4]uint64{ + 0x5a5a5a5a5a5a5a5a, 0x5a5a5a5a5a5a5a5a, 0x5a5a5a5a5a5a5a5a, + 0x5a5a5a5a5a5a5a5a, }, overflow: false, }} @@ -305,35 +357,35 @@ func TestModNScalarBytes(t *testing.T) { in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", expected: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", }, { - name: "group order + 1 (aka 1, overflow in word zero)", + name: "group order + 1 (aka 1, overflow in limb zero)", in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364142", expected: "0000000000000000000000000000000000000000000000000000000000000001", }, { - name: "group order word zero", + name: "group order limb zero", in: "d0364141", expected: "00000000000000000000000000000000000000000000000000000000d0364141", }, { - name: "group order word zero and one", + name: "group order limb zero and one", in: "bfd25e8cd0364141", expected: "000000000000000000000000000000000000000000000000bfd25e8cd0364141", }, { - name: "group order words zero, one, and two", + name: "group order limbs zero, one, and two", in: "af48a03bbfd25e8cd0364141", expected: "0000000000000000000000000000000000000000af48a03bbfd25e8cd0364141", }, { - name: "overflow in word one", + name: "overflow in limb one", in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8dd0364141", expected: "0000000000000000000000000000000000000000000000000000000100000000", }, { - name: "overflow in word two", + name: "overflow in limb two", in: "fffffffffffffffffffffffffffffffebaaedce6af48a03cbfd25e8cd0364141", expected: "0000000000000000000000000000000000000000000000010000000000000000", }, { - name: "overflow in word three", + name: "overflow in limb three", in: "fffffffffffffffffffffffffffffffebaaedce7af48a03bbfd25e8cd0364141", expected: "0000000000000000000000000000000000000001000000000000000000000000", }, { - name: "overflow in word four", + name: "overflow in limb four", in: "ffffffffffffffffffffffffffffffffbaaedce6af48a03bbfd25e8cd0364141", expected: "0000000000000000000000000000000100000000000000000000000000000000", }, { @@ -501,11 +553,11 @@ func TestModNScalarEqualsRandom(t *testing.T) { t.Fatalf("failed equality check\nscalar in: %v", s) } - // Flip a random bit in a random word and ensure it's no longer equal. - randomWord := rng.Int31n(int32(len(s.n))) - randomBit := uint32(1 << uint32(rng.Int31n(32))) + // Flip a random bit in a random limb and ensure it's no longer equal. + randomLimb := rng.Int31n(int32(len(s.n))) + randomBit := uint64(1 << uint64(rng.Int63n(64))) s2 := new(ModNScalar).Set(s) - s2.n[randomWord] ^= randomBit + s2.n[randomLimb] ^= randomBit if s2.Equals(s) { t.Fatalf("failed inequality check\nscalar in: %v", s2) } @@ -531,55 +583,35 @@ func TestModNScalarAdd(t *testing.T) { in2: "0", expected: "1", }, { - name: "group order (aka 0) + 1 (gets reduced, no overflow)", - in1: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141", - in2: "1", - expected: "1", - }, { - name: "group order - 1 + 1 (aka 0, overflow to prime)", + name: "group order-1 + 1 (aka 0, overflow to order)", in1: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", in2: "1", expected: "0", }, { - name: "group order - 1 + 2 (aka 1, overflow in word zero)", - in1: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", - in2: "2", + name: "group order (aka 0) + 1 (aka 1, gets reduced, no overflow)", + in1: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141", + in2: "1", expected: "1", }, { - name: "overflow in word one", - in1: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8bd0364141", - in2: "100000001", + name: "group order-1 + 2 (aka 1, overflow in limb zero)", + in1: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", + in2: "2", expected: "1", }, { - name: "overflow in word two", + name: "overflow in limb one", in1: "fffffffffffffffffffffffffffffffebaaedce6af48a03abfd25e8cd0364141", in2: "10000000000000001", expected: "1", }, { - name: "overflow in word three", - in1: "fffffffffffffffffffffffffffffffebaaedce5af48a03bbfd25e8cd0364141", - in2: "1000000000000000000000001", - expected: "1", - }, { - name: "overflow in word four", + name: "overflow in limb two", in1: "fffffffffffffffffffffffffffffffdbaaedce6af48a03bbfd25e8cd0364141", in2: "100000000000000000000000000000001", expected: "1", }, { - name: "overflow in word five", - in1: "fffffffffffffffffffffffefffffffebaaedce6af48a03bbfd25e8cd0364141", - in2: "10000000000000000000000000000000000000001", - expected: "1", - }, { - name: "overflow in word six", + name: "overflow in limb three", in1: "fffffffffffffffefffffffffffffffebaaedce6af48a03bbfd25e8cd0364141", in2: "1000000000000000000000000000000000000000000000001", expected: "1", - }, { - name: "overflow in word seven", - in1: "fffffffefffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141", - in2: "100000000000000000000000000000000000000000000000000000001", - expected: "1", }, { name: "alternating bits", in1: "a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5", @@ -601,7 +633,7 @@ func TestModNScalarAdd(t *testing.T) { // Ensure the result has the expected value. s1.Add(s2) if !s1.Equals(expected) { - t.Errorf("%s: unexpected result\ngot: %x\nwant: %x", test.name, + t.Errorf("%s: unexpected result\ngot: %v\nwant: %v", test.name, s1, expected) continue } @@ -645,53 +677,6 @@ func TestModNScalarAddRandom(t *testing.T) { } } -// TestAccumulator96Add ensures that the internal 96-bit accumulator used by -// multiplication works as expected for overflow edge cases including overflow. -func TestAccumulator96Add(t *testing.T) { - tests := []struct { - name string // test description - start accumulator96 // starting value of accumulator - in uint64 // value to add to accumulator - expected accumulator96 // expected value of accumulator after addition - }{{ - name: "0 + 0 = 0", - start: accumulator96{[3]uint32{0, 0, 0}}, - in: 0, - expected: accumulator96{[3]uint32{0, 0, 0}}, - }, { - name: "overflow in word zero", - start: accumulator96{[3]uint32{0xffffffff, 0, 0}}, - in: 1, - expected: accumulator96{[3]uint32{0, 1, 0}}, - }, { - name: "overflow in word one", - start: accumulator96{[3]uint32{0, 0xffffffff, 0}}, - in: 0x100000000, - expected: accumulator96{[3]uint32{0, 0, 1}}, - }, { - name: "overflow in words one and two", - start: accumulator96{[3]uint32{0xffffffff, 0xffffffff, 0}}, - in: 1, - expected: accumulator96{[3]uint32{0, 0, 1}}, - }, { - // Start accumulator at 129127208455837319175 which is the result of - // 4294967295 * 4294967295 accumulated seven times. - name: "max result from eight adds of max uint32 multiplications", - start: accumulator96{[3]uint32{7, 4294967282, 6}}, - in: 18446744065119617025, - expected: accumulator96{[3]uint32{8, 4294967280, 7}}, - }} - - for _, test := range tests { - acc := test.start - acc.Add(test.in) - if acc.n != test.expected.n { - t.Errorf("%s: wrong result\ngot: %v\nwant: %v", test.name, acc.n, - test.expected.n) - } - } -} - // TestModNScalarMul ensures that multiplying two scalars together works as // expected for edge cases. func TestModNScalarMul(t *testing.T) { @@ -720,6 +705,16 @@ func TestModNScalarMul(t *testing.T) { in1: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", in2: "2", expected: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd036413f", + }, { + name: "(group order-1) * 3", + in1: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", + in2: "3", + expected: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd036413e", + }, { + name: "(group order-1) * 4", + in1: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", + in2: "4", + expected: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd036413d", }, { name: "(group order-1) * (group order-1)", in1: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", @@ -730,56 +725,76 @@ func TestModNScalarMul(t *testing.T) { in1: "7fffffffffffffffffffffffffffffff5d576e7357a4501ddfe92f46681b20a1", in2: "2", expected: "1", + }, { + name: "highest bit in limb zero squared", + in1: "8000000000000000", + in2: "8000000000000000", + expected: "40000000000000000000000000000000", + }, { + name: "highest bit in limb one squared", + in1: "80000000000000000000000000000000", + in2: "80000000000000000000000000000000", + expected: "4000000000000000000000000000000000000000000000000000000000000000", + }, { + name: "highest bit in limb two squared", + in1: "800000000000000000000000000000000000000000000000", + in2: "800000000000000000000000000000000000000000000000", + expected: "515448c6542dd7f1100b685ccbf26fafc0000000000000000000000000000000", + }, { + name: "highest bit in limb three squared", + in1: "8000000000000000000000000000000000000000000000000000000000000000", + in2: "8000000000000000000000000000000000000000000000000000000000000000", + expected: "2759c7356071a6f179a5fd7916f341f19d0525b0839f3e1e225b3c8519f5f450", + }, { + name: "cross limb product 64x128", + in1: "10000000000000000", + in2: "100000000000000000000000000000000", + expected: "1000000000000000000000000000000000000000000000000", + }, { + name: "cross limb product 64x192", + in1: "10000000000000000", + in2: "1000000000000000000000000000000000000000000000000", + expected: "14551231950b75fc4402da1732fc9bebf", + }, { + name: "cross limb product 128x192", + in1: "100000000000000000000000000000000", + in2: "1000000000000000000000000000000000000000000000000", + expected: "14551231950b75fc4402da1732fc9bebf0000000000000000", }, { name: "group order (aka 0) * 3", in1: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141", in2: "3", expected: "0", }, { - name: "overflow in word eight", + name: "overflow in limb four", in1: "100000000000000000000000000000000", in2: "100000000000000000000000000000000", expected: "14551231950b75fc4402da1732fc9bebf", }, { - name: "overflow in word nine", - in1: "1000000000000000000000000000000000000", - in2: "1000000000000000000000000000000000000", - expected: "14551231950b75fc4402da1732fc9bebf00000000", - }, { - name: "overflow in word ten", + name: "overflow in limb five", in1: "10000000000000000000000000000000000000000", in2: "10000000000000000000000000000000000000000", expected: "14551231950b75fc4402da1732fc9bebf0000000000000000", }, { - name: "overflow in word eleven", - in1: "100000000000000000000000000000000000000000000", - in2: "100000000000000000000000000000000000000000000", - expected: "14551231950b75fc4402da1732fc9bebf000000000000000000000000", - }, { - name: "overflow in word twelve", + name: "overflow in limb six", in1: "1000000000000000000000000000000000000000000000000", in2: "1000000000000000000000000000000000000000000000000", expected: "4551231950b75fc4402da1732fc9bec04551231950b75fc4402da1732fc9bebf", }, { - name: "overflow in word thirteen", - in1: "10000000000000000000000000000000000000000000000000000", - in2: "10000000000000000000000000000000000000000000000000000", - expected: "50b75fc4402da1732fc9bec09d671cd51b343a1b66926b57d2a4c1c61536bda7", - }, { - name: "overflow in word fourteen", + name: "overflow in limb seven", in1: "100000000000000000000000000000000000000000000000000000000", in2: "100000000000000000000000000000000000000000000000000000000", expected: "402da1732fc9bec09d671cd581c69bc59509b0b074ec0aea8f564d667ec7eb3c", }, { - name: "overflow in word fifteen", - in1: "1000000000000000000000000000000000000000000000000000000000000", - in2: "1000000000000000000000000000000000000000000000000000000000000", - expected: "2fc9bec09d671cd581c69bc5e697f5e41f12c33a0a7b6f4e3302b92ea029cecd", + name: "max limb * one", + in1: "ffffffffffffffff", + in2: "1", + expected: "ffffffffffffffff", }, { - name: "double overflow in internal accumulator", - in1: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", - in2: "55555555555555555555555555555554e8e4f44ce51835693ff0ca2ef01215c2", - expected: "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa9d1c9e899ca306ad27fe1945de0242b7f", + name: "max limb * max limb", + in1: "ffffffffffffffff", + in2: "ffffffffffffffff", + expected: "fffffffffffffffe0000000000000001", }, { name: "alternating bits", in1: "a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5", @@ -881,35 +896,35 @@ func TestModNScalarSquare(t *testing.T) { in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364142", expected: "1", }, { - name: "overflow in word eight", + name: "overflow in limb eight", in: "100000000000000000000000000000000", expected: "14551231950b75fc4402da1732fc9bebf", }, { - name: "overflow in word nine", + name: "overflow in limb nine", in: "1000000000000000000000000000000000000", expected: "14551231950b75fc4402da1732fc9bebf00000000", }, { - name: "overflow in word ten", + name: "overflow in limb ten", in: "10000000000000000000000000000000000000000", expected: "14551231950b75fc4402da1732fc9bebf0000000000000000", }, { - name: "overflow in word eleven", + name: "overflow in limb eleven", in: "100000000000000000000000000000000000000000000", expected: "14551231950b75fc4402da1732fc9bebf000000000000000000000000", }, { - name: "overflow in word twelve", + name: "overflow in limb twelve", in: "1000000000000000000000000000000000000000000000000", expected: "4551231950b75fc4402da1732fc9bec04551231950b75fc4402da1732fc9bebf", }, { - name: "overflow in word thirteen", + name: "overflow in limb thirteen", in: "10000000000000000000000000000000000000000000000000000", expected: "50b75fc4402da1732fc9bec09d671cd51b343a1b66926b57d2a4c1c61536bda7", }, { - name: "overflow in word fourteen", + name: "overflow in limb fourteen", in: "100000000000000000000000000000000000000000000000000000000", expected: "402da1732fc9bec09d671cd581c69bc59509b0b074ec0aea8f564d667ec7eb3c", }, { - name: "overflow in word fifteen", + name: "overflow in limb fifteen", in: "1000000000000000000000000000000000000000000000000000000000000", expected: "2fc9bec09d671cd581c69bc5e697f5e41f12c33a0a7b6f4e3302b92ea029cecd", }, { @@ -1015,33 +1030,41 @@ func TestModNScalarNegate(t *testing.T) { in: "1", expected: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", }, { - name: "negation in word one", - in: "0000000000000000000000000000000000000000000000000000000100000000", - expected: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8bd0364141", + name: "two", + in: "2", + expected: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd036413f", + }, { + name: "three", + in: "3", + expected: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd036413e", }, { - name: "negation in word two", + name: "group order - 2", + in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd036413f", + expected: "2", + }, { + name: "group order - 1", + in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", + expected: "1", + }, { + name: "negation in limb one", in: "0000000000000000000000000000000000000000000000010000000000000000", expected: "fffffffffffffffffffffffffffffffebaaedce6af48a03abfd25e8cd0364141", }, { - name: "negation in word three", - in: "0000000000000000000000000000000000000001000000000000000000000000", - expected: "fffffffffffffffffffffffffffffffebaaedce5af48a03bbfd25e8cd0364141", - }, { - name: "negation in word four", + name: "negation in limb two", in: "0000000000000000000000000000000100000000000000000000000000000000", expected: "fffffffffffffffffffffffffffffffdbaaedce6af48a03bbfd25e8cd0364141", }, { - name: "negation in word five", - in: "0000000000000000000000010000000000000000000000000000000000000000", - expected: "fffffffffffffffffffffffefffffffebaaedce6af48a03bbfd25e8cd0364141", - }, { - name: "negation in word six", + name: "negation in limb three", in: "0000000000000001000000000000000000000000000000000000000000000000", expected: "fffffffffffffffefffffffffffffffebaaedce6af48a03bbfd25e8cd0364141", }, { - name: "negation in word seven", - in: "0000000100000000000000000000000000000000000000000000000000000000", - expected: "fffffffefffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141", + name: "half group order (floor)", + in: "7fffffffffffffffffffffffffffffff5d576e7357a4501ddfe92f46681b20a0", + expected: "7fffffffffffffffffffffffffffffff5d576e7357a4501ddfe92f46681b20a1", + }, { + name: "half group order (ceil)", + in: "7fffffffffffffffffffffffffffffff5d576e7357a4501ddfe92f46681b20a1", + expected: "7fffffffffffffffffffffffffffffff5d576e7357a4501ddfe92f46681b20a0", }, { name: "alternating bits", in: "a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5", @@ -1157,33 +1180,29 @@ func TestModNScalarInverseNonConst(t *testing.T) { in: "1", expected: "1", }, { - name: "inverse carry in word one", - in: "0000000000000000000000000000000000000000000000000000000100000000", - expected: "5588b13effffffffffffffffffffffff934e5b00ca8417bf50177f7ba415411a", + name: "two", + in: "2", + expected: "7fffffffffffffffffffffffffffffff5d576e7357a4501ddfe92f46681b20a1", + }, { + name: "group order - 1", + in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", + expected: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", }, { - name: "inverse carry in word two", + name: "group order - 2", + in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd036413f", + expected: "7fffffffffffffffffffffffffffffff5d576e7357a4501ddfe92f46681b20a0", + }, { + name: "inverse carry in limb one", in: "0000000000000000000000000000000000000000000000010000000000000000", expected: "4b0dff665588b13effffffffffffffffa09f710af01555259d4ad302583de6dc", }, { - name: "inverse carry in word three", - in: "0000000000000000000000000000000000000001000000000000000000000000", - expected: "34b9ec244b0dff665588b13effffffffbcff4127932a971a78274c9d74176b38", - }, { - name: "inverse carry in word four", + name: "inverse carry in limb two", in: "0000000000000000000000000000000100000000000000000000000000000000", expected: "50a51ac834b9ec244b0dff665588b13e9984d5b3cf80ef0fd6a23766a3ee9f22", }, { - name: "inverse carry in word five", - in: "0000000000000000000000010000000000000000000000000000000000000000", - expected: "27cfab5e50a51ac834b9ec244b0dff6622f16e85b683d5a059bcd5a3b29d9dff", - }, { - name: "inverse carry in word six", + name: "inverse carry in limb three", in: "0000000000000001000000000000000000000000000000000000000000000000", expected: "897f30c127cfab5e50a51ac834b9ec239c53f268b4700c14f19b9499ac58d8ad", - }, { - name: "inverse carry in word seven", - in: "0000000100000000000000000000000000000000000000000000000000000000", - expected: "6494ef93897f30c127cfab5e50a51ac7b4e8f713e0cddd182234e907286ae6b3", }, { name: "alternating bits", in: "a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5", @@ -1308,21 +1327,33 @@ func TestModNScalarIsOverHalfOrder(t *testing.T) { in: "7fffffffffffffffffffffffffffffff5d576e7357a4501ddfe92f46681b20a1", expected: true, }, { - name: "over half order word one", - in: "7fffffffffffffffffffffffffffffff5d576e7357a4501ddfe92f47681b20a0", + name: "over half order limb one", + in: "7fffffffffffffffffffffffffffffff5d576e7357a4501edfe92f46681b20a0", expected: true, }, { - name: "over half order word two", - in: "7fffffffffffffffffffffffffffffff5d576e7357a4501edfe92f46681b20a0", + name: "group order - 2", + in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd036413f", + expected: true, + }, { + name: "group order - 1", + in: "fffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364140", expected: true, }, { - name: "over half order word three", - in: "7fffffffffffffffffffffffffffffff5d576e7457a4501ddfe92f46681b20a0", + name: "highest limb below half order", + in: "7fffffffffffffffffffffffffffffff00000000000000000000000000000000", + expected: false, + }, { + name: "highest bit set", + in: "8000000000000000000000000000000000000000000000000000000000000000", expected: true, }, { - name: "over half order word seven", - in: "8fffffffffffffffffffffffffffffff5d576e7357a4501ddfe92f46681b20a0", + name: "alternating bits", + in: "a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5", expected: true, + }, { + name: "alternating bits 2", + in: "5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a", + expected: false, }} for _, test := range tests { From 6c0e3cefeec50b52eb5027d80f9e8acca0b0894a Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Thu, 9 Jul 2026 01:28:34 -0500 Subject: [PATCH 11/16] secp256k1: Add scalar reduction proofs. This adds a formal verification proof of the scalar reduction arithmetic in the new 4x64 scalar implementation and updates the proofs README.md accordingly. It includes a formal proof for the following: - 512-bit modular reduction over the group order with saturated 64-bit limbs --- dcrec/secp256k1/internal/proofs/README.md | 9 +- .../proofs/modnscalar_4x64_reduce_prover.py | 232 ++++++++++++++++++ dcrec/secp256k1/modnscalar.go | 4 + 3 files changed, 243 insertions(+), 2 deletions(-) create mode 100644 dcrec/secp256k1/internal/proofs/modnscalar_4x64_reduce_prover.py diff --git a/dcrec/secp256k1/internal/proofs/README.md b/dcrec/secp256k1/internal/proofs/README.md index ba7e919f6..079069a8f 100644 --- a/dcrec/secp256k1/internal/proofs/README.md +++ b/dcrec/secp256k1/internal/proofs/README.md @@ -6,8 +6,9 @@ proofs This consists of formal verification artifacts used to help establish the correctness of the security-critical arithmetic implementations. -The current proofs formally verify properties of the optimized field arithmetic -implementations using [z3](https://github.com/Z3Prover/z3): +The current proofs formally verify properties of the optimized scalar and field +arithmetic implementations using the [Z3 Theorem +Prover](https://github.com/Z3Prover/z3): - Correctness of multi-precision multiplication - Bounds on intermediate values @@ -45,3 +46,7 @@ optimized implementations preserve the required arithmetic invariants. - [Field Modular Reduction](field_4x64_reduce_prover.py) Provides formal verification of the reduction of a 512-bit value represented using saturated 64-bit limbs modulo the secp256k1 prime. + +- [Scalar Modular Reduction](modnscalar_4x64_reduce_prover.py) + Provides formal verification of the reduction of a 512-bit value represented + using saturated 64-bit limbs modulo the secp256k1 group order. diff --git a/dcrec/secp256k1/internal/proofs/modnscalar_4x64_reduce_prover.py b/dcrec/secp256k1/internal/proofs/modnscalar_4x64_reduce_prover.py new file mode 100644 index 000000000..d486c55ef --- /dev/null +++ b/dcrec/secp256k1/internal/proofs/modnscalar_4x64_reduce_prover.py @@ -0,0 +1,232 @@ +# Copyright (c) 2026 The Decred developers +# Copyright (c) 2026 Dave Collins +# Use of this source code is governed by an ISC +# license that can be found in the LICENSE file. + +from z3_proof_helpers import * + +# ------- +# Inputs. +# ------- + +x = BitVecs('x0 x1 x2 x3 x4 x5 x6 x7', 64) + +# NOTE: The code in this section reuses the variable names from the Go code +# being proven for easy cross referencing. However, the comments refer to them +# by shorter mathematical variables for clarity. + +# Define N. +orderLimb0 = BitVecVal(0xbfd25e8cd0364141, 64) +orderLimb1 = BitVecVal(0xbaaedce6af48a03b, 64) +orderLimb2 = BitVecVal(0xfffffffffffffffe, 64) +orderLimb3 = BitVecVal(0xffffffffffffffff, 64) +order = Concat(orderLimb3, orderLimb2, orderLimb1, orderLimb0) + +# Define 2N. +twiceOrderLimb0 = BitVecVal(0x7fa4bd19a06c8282, 64) +twiceOrderLimb1 = BitVecVal(0x755db9cd5e914077, 64) +twiceOrderLimb2 = BitVecVal(0xfffffffffffffffd, 64) +twiceOrderLimb3 = BitVecVal(0xffffffffffffffff, 64) +twiceOrderLimb4 = ONE +twiceOrder = Concat( + twiceOrderLimb4, twiceOrderLimb3, twiceOrderLimb2, + twiceOrderLimb1, twiceOrderLimb0 +) + +# Define 2**256 - N. +orderComplementLimb0 = BitVecVal(0x402da1732fc9bebf, 64) +orderComplementLimb1 = BitVecVal(0x4551231950b75fc4, 64) +orderComplementLimb2 = ONE +orderComplement = Concat(orderComplementLimb2, orderComplementLimb1, orderComplementLimb0) + +# -------------------------------------------------------------------------- +# Prove self-consistency of the constants before trusting anything that uses +# them. +# -------------------------------------------------------------------------- + +# Prove twiceOrder = 2N. +order320 = ZeroExt(64, order) +prove(twiceOrder == order320 << 1, "twiceOrder != 2N") + +# -------------------------------------------------------------------------- +# Lemma R-identity: N + c == 2**256 +# +# This lemma justifies replacing multiples of 2**256 with multiples of c when +# working modulo N: +# +# N + c ≡ 2**256 (mod N) +# => 0 + c ≡ 2**256 (mod N) +# => c ≡ 2**256 (mod N) +# -------------------------------------------------------------------------- +orderComplement320 = ZeroExt(128, orderComplement) +prove(order320+orderComplement320 == BitVecVal(1, 320)<<256, + "Lemma R-identity: N+c != 2**256") + +# --------------- +# Model the code. +# --------------- + +discards = [] + +# first reduction: 512 bits -> 385 bits. +h0, t0 = mul64(x[4], orderComplementLimb0) +h1, t1 = mul64(x[4], orderComplementLimb1) +t1, carry = add64(t1, h0, ZERO) +t2, carry = add64(x[4], h1, carry) +t3 = carry + +h0, l0 = mul64(x[5], orderComplementLimb0) +h1, l1 = mul64(x[5], orderComplementLimb1) +t1, carry = add64(t1, l0, ZERO) +t2, carry = add64(t2, h0, carry) +t3, discarded = add64(t3, h1, carry) +discards.append(discarded) +t2, carry = add64(t2, l1, ZERO) +t3, carry = add64(t3, x[5], carry) +t4 = carry + +h0, l0 = mul64(x[6], orderComplementLimb0) +h1, l1 = mul64(x[6], orderComplementLimb1) +t2, carry = add64(t2, l0, ZERO) +t3, carry = add64(t3, h0, carry) +t4, discarded = add64(t4, h1, carry) +discards.append(discarded) +t3, carry = add64(t3, l1, ZERO) +t4, carry = add64(t4, x[6], carry) +t5 = carry + +h0, l0 = mul64(x[7], orderComplementLimb0) +h1, l1 = mul64(x[7], orderComplementLimb1) +t3, carry = add64(t3, l0, ZERO) +t4, carry = add64(t4, h0, carry) +t5, discarded = add64(t5, h1, carry) +discards.append(discarded) +t4, carry = add64(t4, l1, ZERO) +t5, carry = add64(t5, x[7], carry) +t6 = carry + +t0, carry = add64(t0, x[0], ZERO) +t1, carry = add64(t1, x[1], carry) +t2, carry = add64(t2, x[2], carry) +t3, carry = add64(t3, x[3], carry) +t4, carry = add64(t4, ZERO, carry) +t5, carry = add64(t5, ZERO, carry) +t6, discarded = add64(t6, ZERO, carry) +discards.append(discarded) + +# Save for analysis. +t6_saved_1 = t6 +t_full_redux_1 = Concat(ZERO, t6, t5, t4, t3, t2, t1, t0) +first_redux_input_high = Concat(x[7], x[6], x[5], x[4]) + +# second reduction: 385 bits -> 258 bits. +h0, l0 = mul64(t4, orderComplementLimb0) +h1, l1 = mul64(t4, orderComplementLimb1) +l2 = t4 +t0, carry = add64(t0, l0, ZERO) +t1, carry = add64(t1, h0, carry) +t2, carry = add64(t2, h1, carry) +t3, carry = add64(t3, ZERO, carry) +t4 = carry +t1, carry = add64(t1, l1, ZERO) +t2, carry = add64(t2, l2, carry) +t3, carry = add64(t3, ZERO, carry) +t4, discarded = add64(t4, ZERO, carry) +discards.append(discarded) + +h0, l0 = mul64(t5, orderComplementLimb0) +h1, l1 = mul64(t5, orderComplementLimb1) +t1, carry = add64(t1, l0, ZERO) +t2, carry = add64(t2, h0, carry) +t3, carry = add64(t3, h1, carry) +t4, discarded = add64(t4, ZERO, carry) +discards.append(discarded) +t2, carry = add64(t2, l1, ZERO) +t3, carry = add64(t3, t5, carry) +t4, discarded = add64(t4, ZERO, carry) +discards.append(discarded) + +t2, carry = add64(t2, t6*orderComplementLimb0, ZERO) +t3, carry = add64(t3, t6*orderComplementLimb1, carry) +t4, discarded = add64(t4, t6, carry) +discards.append(discarded) + +# Save intermediate results for analysis. +t3_carry_2 = carry +t3_saved_2 = t3 +t4_saved_2 = t4 + +# third reduction: max 258 bits -> t < 2N. +t0, carry = add64(t0, t4*orderComplementLimb0, ZERO) +t1, carry = add64(t1, t4*orderComplementLimb1, carry) +t2, carry = add64(t2, t4, carry) +t3, carry = add64(t3, ZERO, carry) +t4 = carry + +# final reduction: t < 2N -> t < N +s0, borrow = sub64(t0, orderLimb0, ZERO) +s1, borrow = sub64(t1, orderLimb1, borrow) +s2, borrow = sub64(t2, orderLimb2, borrow) +s3, borrow = sub64(t3, orderLimb3, borrow) +_, borrow = sub64(t4, ZERO, borrow) +r0 = If(borrow == ONE, t0, s0) +r1 = If(borrow == ONE, t1, s1) +r2 = If(borrow == ONE, t2, s2) +r3 = If(borrow == ONE, t3, s3) + +# ------ +# Proofs. +# ------- + +# Discarded carries are never set. +prove_no_discarded_carries(discards) + +# Top limb after first reduction is max 1 (t6 ≤ 1). +# +# This along with the proof there is no discarded carry on t6 after the first +# reduction proves the value does not exceed 385 bits (as expected since the +# complement is 129 bits and there were a max of 256 bits in the upper limbs +# before the first reduction, so 256 + 129 = 385). +prove(ULE(t6_saved_1, ONE), "top limb after 1st reduction > 1") + +# Top limb after second reduction is never exceeds 3 (t4 ≤ 3). +# +# This along with the proof there is no discarded carry on t4 after the second +# reduction proves the value does not exceed 258 bits (as expected because the +# complement is 129 bits and there were 129 max bits left in the upper limbs +# after the first reduction, so 129 + 129 = 258). +prove(ULE(t4_saved_2, BitVecVal(3, 64)), "top limb after 2nd reduction > 3") + +# Top limb after second reduction times complement limb 0 never exceeds uint64. +# ((t4*c0)>>64 == 0) +# +# This proves the regular 64-bit product used in the third reduction does not +# overflow. +t4c0hi, _ = mul64(t4_saved_2, orderComplementLimb0) +prove(t4c0hi == ZERO, + "high word of top limb from 2nd reduction times complement limb 0 != 0") + +# Top limb after second reduction times complement limb 1 never exceeds uint64. +# ((t4*c1)>>64 == 0) +# +# This proves the regular 64-bit product used in the third reduction does not +# overflow. +t4c1hi, _ = mul64(t4_saved_2, orderComplementLimb1) +prove(t4c1hi == ZERO, + "high word of top limb from 2nd reduction times complement limb 1 != 0") + +# Top limb after third reduction is max 1 (t4 ≤ 1). +# +# This proof is not strictly necessary since the bound is actually tigher per +# the next proof which implicitly proves this fact. However, it is fast to +# prove and including it provides nice symmetry with the previous proofs for the +# bounds of the top limb after each reduction. +prove(ULE(t4, ONE), "top limb after 3rd reduction > 1") + +# Value after third reduction is less than twice the group order (t < 2N). +t = Concat(t4, t3, t2, t1, t0) +prove(ULT(t, twiceOrder), "value after 3rd reduction >= 2N") + +# Fully reduced result is less than the group order (r < N). +r = Concat(r3, r2, r1, r0) +prove(ULT(r, order), "fully reduced result >= N") diff --git a/dcrec/secp256k1/modnscalar.go b/dcrec/secp256k1/modnscalar.go index 6ccaf6306..55356dad1 100644 --- a/dcrec/secp256k1/modnscalar.go +++ b/dcrec/secp256k1/modnscalar.go @@ -387,6 +387,10 @@ func (s *ModNScalar) Add(val *ModNScalar) *ModNScalar { // scalar64Reduce512 reduces a 512-bit little-endian limb array modulo the group // order in constant time and stores the result in r. func scalar64Reduce512(r *[4]uint64, x *[8]uint64) { + // The intermediate bounds and carry assumptions used by this algorithm have + // been formally verified. The verification artifacts are available in + // internal/proofs. + // The overall strategy employed here is: // 1) Start with the full unreduced 512-bit product of the two scalars. // 2) Reduce the result modulo the group order via Crandall reduction From dbce45af1ef0ac21a04fdaaa9f039b669a6a396c Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Sun, 19 Jul 2026 08:06:40 -0500 Subject: [PATCH 12/16] secp256k1: Add width-5 windowed NAF recoder. With the ultimate goal of optimizing scalar multiplication by reducing the average number of point additions, this commit adds code to convert a balanced scalar representative into width-5 windowed non-adjacent form (wNAF). The current scalar multiplication implementation uses ordinary non-adjacent form (NAF), which corresponds to a window width of 2. This commit introduces only the width-5 wNAF recoder and does not yet integrate it into scalar multiplication. Future commits will add comprehensive tests and integrate the new implementation. --- dcrec/secp256k1/curve.go | 228 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 228 insertions(+) diff --git a/dcrec/secp256k1/curve.go b/dcrec/secp256k1/curve.go index 750b57a32..fe93eb8da 100644 --- a/dcrec/secp256k1/curve.go +++ b/dcrec/secp256k1/curve.go @@ -943,6 +943,234 @@ func splitK(k *ModNScalar) (ModNScalar, ModNScalar) { return k1, k2 } +// wNAFWidth is the window width of the NAF representation used in scalar +// multiplication. +// +// It is important to note that this constant primarily exists to improve code +// readability, to parameterize parts of the scalar multiplication +// implementation, and to allow tests to assert invariants. +// +// However, the overall implementation is intentionally not fully parameterized +// based on this constant because it is carefully optimized for this specific +// window width and those optimizations involve taking advantage of known bounds +// that an arbitrary window size would violate. +// +// Choosing a new window size involves a variety of tradeoffs that must be +// carefully analyzed. For example, recoding costs, precomputation costs, and +// average density. +const wNAFWidth = 5 + +// wnaf5RecodeCodes maps the low five bits of an integer to the encoded +// representative of the width-5 wNAF digit that should be emitted. +// +// These values are not the signed digits themselves. Rather, each is an +// encoded representative that serves directly as an index into precomputed +// tables. See [wnafScalar] for the encoding details. +// +// For each odd residue, the represented signed digit is the unique value in +// {±1, ±3, ±5, ±7, ±9, ±11, ±13, ±15} that is congruent to the residue modulo +// 32: +// +// residue ≡ digit (mod 32) +// +// Thus, subtracting the represented signed digit always produces an integer +// that is divisible by 32. This allows the optimized recoding algorithm to +// skip five one-bit iterations of the canonical algorithm at once. +// +// Even residues always map to zero because they are already divisible by two +// and therefore correspond only to implicit zero digits. +// +// The corresponding arithmetic adjustment for each encoded representative is +// provided by [wnaf5RecodeAdjustments]. +var wnaf5RecodeCodes = [32]uint8{ + 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, 0, 8, + 0, 16, 0, 15, 0, 14, 0, 13, 0, 12, 0, 11, 0, 10, 0, 9, +} + +// wnaf5RecodeAdjustments maps the low five bits of an integer to the signed +// adjustment that must be added to the working integer before shifting during +// the width-5 wNAF recoding algorithm. +// +// While [wnaf5RecodeCodes] produces the encoded representative that is stored +// in the resulting wNAF encoding, this table provides the corresponding +// arithmetic adjustment. The represented signed digit d satisfies: +// +// residue ≡ d (mod 32) +// +// and the recoding algorithm updates the working integer via: +// +// n = (n - d) / 2 +// +// Rather than subtracting d directly, this table stores -d encoded as an +// unsigned 64-bit two's complement value. Negative adjustments therefore +// appear as values such as: +// +// -1 -> ^uint64(0) +// -3 -> ^uint64(2) +// -5 -> ^uint64(4) +// -7 -> ^uint64(6) +// -9 -> ^uint64(8) +// -11 -> ^uint64(10) +// -13 -> ^uint64(12) +// -15 -> ^uint64(14) +// +// The adjustment is sign-extended across the upper limbs during the multiword +// addition which allows the hot loop to remain branch free. +var wnaf5RecodeAdjustments = [32]uint64{ + 0, ^uint64(0), 0, ^uint64(2), 0, ^uint64(4), 0, ^uint64(6), + 0, ^uint64(8), 0, ^uint64(10), 0, ^uint64(12), 0, ^uint64(14), + 0, 15, 0, 13, 0, 11, 0, 9, 0, 7, 0, 5, 0, 3, 0, 1, +} + +// wnafScalar represents a non-negative integer less than 2^129 encoded in +// width-5 windowed non-adjacent form (wNAF). +// +// wNAF is a signed-digit representation where the valid digits are 0, ±1, ±3, +// ±5, ±7, ±9, ±11, ±13, and ±15. +// +// Each entry corresponds to one bit position and encodes the signed digit. +// +// The encoding from code to digit is: +// 0 -> 0 +// 1 -> +1 +// 2 -> +3 +// 3 -> +5 +// 4 -> +7 +// 5 -> +9 +// 6 -> +11 +// 7 -> +13 +// 8 -> +15 +// 9 -> -1 +// 10 -> -3 +// 11 -> -5 +// 12 -> -7 +// 13 -> -9 +// 14 -> -11 +// 15 -> -13 +// 16 -> -15 +// +// The encoding is intentionally not using a more compact form so that it serves +// directly as an index into precomputed tables which allows simplification of +// hot paths. +// +// Formally, letting c denote the stored code value, the decoded signed digit +// d(c) is given by the piecewise function: +// +// {0, where c = 0 +// d(c) = {2c - 1, where 0 < c ≤ 8 +// {17 - 2c, where 8 < c ≤ 16 +// +// The array contains one extra entry because a recoding may produce a carry +// into an additional most-significant digit. +type wnafScalar struct { + codes [130]uint8 + bits uint8 +} + +// wnaf takes a non-negative integer less than 2^129 and returns its width-5 +// windowed non-adjacent form (wNAF) which is a unique signed-digit +// representation such that non-zero digits are separated by at least 4 zeroes. +// See [wnafScalar] for details on how the representation is encoded and how to +// interpret it. +// +// Width-5 wNAF is useful because, on average, only about one in every six +// digits is non-zero. +// +// This property is particularly beneficial for optimizing elliptic curve point +// multiplication because it greatly reduces the number of point additions at +// the cost of precomputing a few odd multiples of the point and, in the worst +// case, one additional point doubling due to a carry introduced by the +// recoding. This is an excellent tradeoff because subtraction of points has +// the same computational complexity as addition of points and point doubling is +// faster than both. +func wnaf(k *ModNScalar) wnafScalar { + const ( + // This is intentionally not using the package constant [wNAFWidth] for + // the window size for the reasons stated by its documentation. + // + // In particular, this implementation is carefully optimized for this + // specific width and involves assumptions that an arbitrary window size + // might violate. + // + // Using a separate constant helps make it clear that updating the + // window size here requires extra care to assert correctness. + windowSize = 5 + windowMask = (1 << windowSize) - 1 + ) + + // The arithmetic for wNAF recoding is performed over the ordinary integers, + // not modulo the curve order. Moreover, this method is required to be + // called with the scalar's balanced integer representative, which is known + // to fit within 129 bits. + // + // Consequently, the scalar is first converted to a uint192 with three + // 64-bit limbs where the top limb is at most 1. + t0, t1, t2 := k.n[0], k.n[1], k.n[2] + + // This implementation is based on the standard width-w NAF recoding + // algorithm presented as Algorithm 3.35 in [GECC]. However, it has been + // modified to skip runs of zero digits instead of processing one bit at a + // time. + // + // The optimizations exploit the fact that runs of zero digits are implicit + // in the output representation. + // + // Also, note that the last non-zero bit is initialized to the maximum value + // so that adding 1 at the end to account for the exclusive endpoint wraps + // around to 0 when no digits are emitted. + var result wnafScalar + var c uint64 + var bit uint8 + var lastNonZeroBit = ^uint8(0) + for t0|t1|t2 != 0 { + // The next five iterations of the canonical bit-by-bit algorithm would + // emit only zero when the low window is zero, so skip directly to the + // next window boundary in that case. + if residue := uint8(t0 & windowMask); residue != 0 { + // Skip the zero digits that the canonical bit-by-bit algorithm + // would emit before reaching the next odd value. + // + // Since 0 < residue < 32, shift is guaranteed to satisfy the + // following bounds: + // + // 0 ≤ shift < 5 + shift := uint8(bits.TrailingZeros8(residue)) + t0 = t0>>shift | t1<<(64-shift) + t1 = t1>>shift | t2<<(64-shift) + t2 >>= shift + bit += shift + + residue = uint8(t0 & windowMask) + result.codes[bit] = wnaf5RecodeCodes[residue] + lastNonZeroBit = bit + + // The selected digit is congruent to the low five bits modulo 32. + // Therefore, adding the stored adjustment (which represents the + // negation of that digit) makes the value divisible by 32. + // + // Negative adjustments are stored in two's complement. Sign + // extending the value across the upper limbs allows the multiword + // addition to operate without branches. + adjustment := wnaf5RecodeAdjustments[residue] + signExtended := uint64(int64(adjustment) >> 63) + t0, c = bits.Add64(t0, adjustment, 0) + t1, c = bits.Add64(t1, signExtended, c) + t2, _ = bits.Add64(t2, signExtended, c) + } + + // Divide by 32 to prepare for the next window. + // + // The adjustment above guarantees the current value is divisible by 32. + t0 = t0>>windowSize | t1<<(64-windowSize) + t1 = t1>>windowSize | t2<<(64-windowSize) + t2 >>= windowSize + bit += windowSize + } + + result.bits = lastNonZeroBit + 1 + return result +} + // nafScalar represents a positive integer up to a maximum value of 2^256 - 1 // encoded in non-adjacent form. // From 5d91d3443f30dbe671345b1209fb257864545fda Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Sun, 19 Jul 2026 08:06:42 -0500 Subject: [PATCH 13/16] secp256k1: Add wNAF tests. This adds comprehensive tests for the new wNAF implementation to help verify correctness. The tests assert the resulting encoding satisfies the required mathematical properties and reconstructs the original value for both deterministic edge cases and randomized balanced scalar representatives. The deterministic edge cases exercise enough small values to exhaustively test every residue modulo 2^w multiple times, and include additional edge cases such as carries introduced by the recoding, maximum values, and alternating bit patterns. --- dcrec/secp256k1/curve_test.go | 175 ++++++++++++++++++++++++++++++++++ 1 file changed, 175 insertions(+) diff --git a/dcrec/secp256k1/curve_test.go b/dcrec/secp256k1/curve_test.go index d3eb497f2..607e1ce32 100644 --- a/dcrec/secp256k1/curve_test.go +++ b/dcrec/secp256k1/curve_test.go @@ -753,6 +753,181 @@ func TestNAFRandom(t *testing.T) { } } +// decodeWNAFDigit converts the provided wNAF code to the signed digit it +// represents. +func decodeWNAFDigit(code uint8) int8 { + const ( + positiveCodeLimit = 1 << (wNAFWidth - 2) + maxCode = 1 << (wNAFWidth - 1) + ) + switch { + case code == 0: + return 0 + case code <= positiveCodeLimit: + return 2*int8(code) - 1 + } + return (maxCode + 1) - 2*int8(code) +} + +// checkWNAFEncoding returns an error if the provided windowed NAF encoding does +// not adhere to the encoding requirements or does not reconstruct the provided +// original value. +func checkWNAFEncoding(s *wnafScalar, origValue *big.Int) error { + // Ensure the reported number of used bits does not exceed the size of the + // array. + if int(s.bits) > len(s.codes) { + return fmt.Errorf("bits %d > %d", s.bits, len(s.codes)) + } + + // wNAF must not have a leading zero and the individual codes must not + // exceed the max allowed code for the window width. + codes := s.codes[:s.bits] + if len(codes) > 0 && codes[len(codes)-1] == 0 { + return fmt.Errorf("leading zero in encoding: %v", codes) + } + const maxCode = 1 << (wNAFWidth - 1) + for _, code := range codes { + if code > maxCode { + return fmt.Errorf("found code %d > %d", code, maxCode) + } + } + + // Ensure each non-zero digit is separated by at least the width of the + // window. + const minDistance = wNAFWidth + prevNonZeroBit := -minDistance + for bit, code := range codes { + if code == 0 { + continue + } + + if distance := bit - prevNonZeroBit; distance < minDistance { + return fmt.Errorf("non-zero digits at bit pos %d and %d are only "+ + "%d bits apart", prevNonZeroBit, bit, distance) + } + + prevNonZeroBit = bit + } + + // Reconstruct the represented value and ensure it matches the original + // value. + sum := new(big.Int) + for bit, code := range codes { + digit := decodeWNAFDigit(code) + term := big.NewInt(int64(digit)) + term.Lsh(term, uint(bit)) + sum.Add(sum, term) + } + if origValue.Cmp(sum) != 0 { + return fmt.Errorf("failed to reconstruct orig value: got %v, want %v", + sum, origValue) + } + + return nil +} + +// TestWNAF ensures encoding various edge cases and values to width-w windowed +// non-adjacent form, where w is [wNAFWidth], produces valid results. +func TestWNAF(t *testing.T) { + type wnafTest struct { + name string // test description + in string // hex encoded test value + } + extraTests := []wnafTest{{ + name: "leading zeroes", + in: "002f20569b90697ad471c1be6107814f", + }, { + name: "highest allowed bit only", + in: "100000000000000000000000000000000", + }, { + name: "largest 129-bit value", + in: "1ffffffffffffffffffffffffffffffff", + }, { + name: "130 bits when NAF encoded", + in: "1f0000000000000000000000000000001", + }, { + name: "alternating bits (0xa5)", + in: "a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5", + }, { + name: "alternating bits (0x5a)", + in: "5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a5a", + }, { + name: "all ones", + in: "ffffffffffffffffffffffffffffffff", + }, { + name: "first term of balanced length-two representation #1", + in: "b776e53fb55f6b006a270d42d64ec2b1", + }, { + name: "second term balanced length-two representation #1", + in: "d6cc32c857f1174b604eefc544f0c7f7", + }, { + name: "first term of balanced length-two representation #2", + in: "45c53aa1bb56fcd68c011e2dad6758e4", + }, { + name: "second term of balanced length-two representation #2", + in: "a2e79d200f27f2360fba57619936159b", + }} + + // Add tests for the first small values since they collectively exercise + // every possible residue modulo 2^w, where w is the window width, multiple + // times. Then append the extra deterministic tests for edge cases. + const smallValues = 2*(1< Date: Sun, 19 Jul 2026 08:06:43 -0500 Subject: [PATCH 14/16] secp256k1: Add wNAF benchmark. --- dcrec/secp256k1/curve_bench_test.go | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/dcrec/secp256k1/curve_bench_test.go b/dcrec/secp256k1/curve_bench_test.go index 5100ed896..5bcd50e17 100644 --- a/dcrec/secp256k1/curve_bench_test.go +++ b/dcrec/secp256k1/curve_bench_test.go @@ -229,6 +229,18 @@ func BenchmarkNAF(b *testing.B) { } } +// BenchmarkWNAF benchmarks conversion of a non-negative integer into its +// width-w windowed non-adjacent form representation, where w is [wNAFWidth]. +func BenchmarkWNAF(b *testing.B) { + k := mustModNScalar("8447e288d34b360bc885cb8ce7c00575") + + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + wnaf(k) + } +} + // BenchmarkJacobianPointEquivalency benchmarks determining if two Jacobian // points represent the same affine point. func BenchmarkJacobianPointEquivalency(b *testing.B) { From 4ee1c16ed28c61b5925b3f034ea94393a6854550 Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Sun, 19 Jul 2026 08:06:45 -0500 Subject: [PATCH 15/16] secp256k1: Integrate width-5 wNAF in scalar mult. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This updates the non-constant-time scalar multiplication implementation to use the new width-5 wNAF recoding. In particular, because wNAF requires precomputed odd multiples that were not needed by the previous implementation, this introduces a new dedicated precomputed table type along with a new helper to construct the required multiples. The previous implementation only needed to swap between two precomputed points depending on the scalar signs. However, since the new wNAF implementation uses larger precomputed tables, it instead tracks whether each decomposed scalar was negated and applies the corresponding sign adjustment to the appropriate half of each table. The following benchmark shows a before and after comparison of scalar multiplication: name old time/op new time/op delta --------------------------------------------------------------------------- ScalarMultNonConst 94.4µs ± 1% 84.3µs ± 1% -10.75% (p=0.000 n=10+10) --- dcrec/secp256k1/curve.go | 243 ++++++++++++++++++++++++++------------- 1 file changed, 163 insertions(+), 80 deletions(-) diff --git a/dcrec/secp256k1/curve.go b/dcrec/secp256k1/curve.go index fe93eb8da..9c2830264 100644 --- a/dcrec/secp256k1/curve.go +++ b/dcrec/secp256k1/curve.go @@ -1307,6 +1307,111 @@ func naf(k []byte) nafScalar { return result } +// wNAFPrecompTableSize is the size of the wNAF precomputed point table. +const wNAFPrecompTableSize = 1 << (wNAFWidth - 1) + +// wNAFPrecompTable houses the precomputed odd point multiples used during wNAF +// scalar multiplication. +// +// The first half of the table (starting at index 1) contains the positive odd +// multiples: +// +// P, 3P, 5P, ... +// +// While the second half of the table contains the corresponding negatives: +// +// -P, -3P, -5P, ... +// +// Entry 0 is intentionally left unused so the encoded wNAF digit can be used +// directly as an array index. +type wNAFPrecompTable [wNAFPrecompTableSize + 1]JacobianPoint + +// wNAFNegateStart returns the index of the half of a precomputed table whose +// points must be negated in order to match the sign adjustment made to the +// corresponding scalar. +// +// When the scalar has been negated, the first half is negated so it represents +// the negative odd multiples while the second half remains positive. +// Otherwise, the opposite arrangement is used. +func wNAFNegateStart(isScalarNegated bool) int { + if isScalarNegated { + return 1 + } + return wNAFPrecompTableSize/2 + 1 +} + +// buildOddPrecomputeTables constructs the odd multiple tables for both P and +// φ(P) and arranges the positive and negative representatives so they can be +// indexed directly by encoded wNAF digits. +func buildOddPrecomputeTables(p *JacobianPoint, pTable, phiTable *wNAFPrecompTable, k1Neg, k2Neg bool) { + // oddMultiplesPerHalf is the total number of odd multiples in each half of + // the table. + const oddMultiplesPerHalf = wNAFPrecompTableSize / 2 + + // Build the positive odd multiples of P. + // + // Width-w wNAF only ever emits odd signed digits, so only odd multiples of + // the point need to be precomputed. + // + // Both positive and negative entries are ultimately needed. The second + // half is initially a copy of the positive half and the appropriate half is + // negated below depending on the sign of k1. + // + // The negation is deferred until after the φ(P) table is built because that + // table is constructed by applying φ to the corresponding positive odd + // multiples of P. + var twoP JacobianPoint + pTable[1].Set(p) + DoubleNonConst(&pTable[1], &twoP) + for i := 1; i < oddMultiplesPerHalf; i++ { + AddNonConst(&twoP, &pTable[i], &pTable[i+1]) + } + copy(pTable[oddMultiplesPerHalf+1:], pTable[1:oddMultiplesPerHalf+1]) + + // Build the corresponding odd multiples of φ(P). + // + // As above, both positive and negative entries are ultimately needed. The + // sign of this table is chosen independently from the P table because it + // depends only on the sign of k2. + // + // Note that φ is a group homomorphism, so: + // + // φ(n*P) = n*φ(P) + // + // This allows the second table to be obtained by applying φ to each entry + // of the first table instead of computing each odd multiple independently. + // + // NOTE: φ(x,y) = (βx,y). The Jacobian z coordinates are the same, so this + // math goes through. + for i := 1; i <= oddMultiplesPerHalf; i++ { + phiTable[i].Set(&pTable[i]) + phiTable[i].X.Mul(endoBeta) + } + copy(phiTable[oddMultiplesPerHalf+1:], phiTable[1:oddMultiplesPerHalf+1]) + + // Negate whichever half of the table corresponds to the sign adjustment + // that was applied to k1. + // + // This maintains the invariant: + // + // k1*P = -k1*-P + negateStart := wNAFNegateStart(k1Neg) + for i := negateStart; i < negateStart+oddMultiplesPerHalf; i++ { + pTable[i].Y.Negate(1).Normalize() + } + + // Apply the same transformation independently to the φ(P) table based on + // the sign adjustment that was applied to k2. + // + // Similarly, this maintains the invariant: + // + // k2*φ(P) = -k2*-φ(P) + negateStart = wNAFNegateStart(k2Neg) + for i := negateStart; i < negateStart+oddMultiplesPerHalf; i++ { + phiTable[i].Y.Negate(1).Normalize() + } +} + // ScalarMultNonConst multiplies k*P where k is a scalar modulo the curve order // and P is a point in Jacobian projective coordinates and stores the result in // the provided Jacobian point. @@ -1354,25 +1459,6 @@ func ScalarMultNonConst(k *ModNScalar, point, result *JacobianPoint) { // See section 3.5 in [GECC] for a more rigorous treatment. // ------------------------------------------------------------------------- - // Per above, the main equation here to remember is: - // k*P = k1*P + k2*φ(P) - // - // p1 below is P in the equation while p2 is φ(P) in the equation. - // - // NOTE: φ(x,y) = (β*x,y). The Jacobian z coordinates are the same, so this - // math goes through. - // - // Also, calculate -p1 and -p2 for use in the NAF optimization. - p1, p1Neg := new(JacobianPoint), new(JacobianPoint) - p1.Set(point) - p1Neg.Set(p1) - p1Neg.Y.Negate(1).Normalize() - p2, p2Neg := new(JacobianPoint), new(JacobianPoint) - p2.Set(p1) - p2.X.Mul(endoBeta).Normalize() - p2Neg.Set(p2) - p2Neg.Y.Negate(1).Normalize() - // Decompose k into k1 and k2 such that k = k1 + k2*λ (mod n) where k1 and // k2 are around half the bit length of k in order to halve the number of EC // operations. @@ -1384,94 +1470,91 @@ func ScalarMultNonConst(k *ModNScalar, point, result *JacobianPoint) { // which means that when they would otherwise be a small negative magnitude // they will instead be a large positive magnitude. Since the goal is for // the scalars to have a small magnitude to achieve a performance boost, use - // their negation when they are greater than the half order of the group and - // flip the positive and negative values of the corresponding point that - // will be multiplied by to compensate. + // their negation when they are greater than the half order of the group. + // + // In order to compensate, the positive and negative values of the + // corresponding points that will be multiplied by are flipped later. // // In other words, transform the calc when k1 is over the half order to: // k1*P = -k1*-P // // Similarly, transform the calc when k2 is over the half order to: // k2*φ(P) = -k2*-φ(P) + var k1Neg, k2Neg bool k1, k2 := splitK(k) if k1.IsOverHalfOrder() { k1.Negate() - p1, p1Neg = p1Neg, p1 + k1Neg = true } if k2.IsOverHalfOrder() { k2.Negate() - p2, p2Neg = p2Neg, p2 + k2Neg = true } - // Convert k1 and k2 into their NAF representations since NAF has a lot more - // zeros overall on average which minimizes the number of required point - // additions in exchange for a mix of fewer point additions and subtractions - // at the cost of one additional point doubling. + // Per above, the main equation here to remember is: + // k*P = k1*P + k2*φ(P) + // + // The simultaneous multiplication therefore needs two sets of precomputed + // odd multiples: + // + // P, 3P, 5P, ... + // + // and: + // + // φ(P), 3φ(P), 5φ(P), ... + var pPrecomps, phiPrecomps wNAFPrecompTable + buildOddPrecomputeTables(point, &pPrecomps, &phiPrecomps, k1Neg, k2Neg) + + // Convert k1 and k2 into their windowed NAF representations since they have + // a lot more zeroes overall on average which greatly reduces the number of + // point additions at the cost of precomputing a few odd multiples of the + // point and, in the worst case, one additional point doubling due to a + // carry introduced by the recoding. // // This is an excellent tradeoff because subtraction of points has the same // computational complexity as addition of points and point doubling is // faster than both. // // Concretely, on average, 1/2 of all bits will be non-zero with the normal - // binary representation whereas only 1/3rd of the bits will be non-zero - // with NAF. - // - // The Pos version of the bytes contain the +1s and the Neg versions contain - // the -1s. - k1Bytes, k2Bytes := k1.Bytes(), k2.Bytes() - k1NAF, k2NAF := naf(k1Bytes[:]), naf(k2Bytes[:]) - k1PosNAF, k1NegNAF := k1NAF.Pos(), k1NAF.Neg() - k2PosNAF, k2NegNAF := k2NAF.Pos(), k2NAF.Neg() - k1Len, k2Len := len(k1PosNAF), len(k2PosNAF) + // binary representation whereas only 1/6 of the bits will be non-zero with + // width-5 wNAF. + k1NAF, k2NAF := wnaf(&k1), wnaf(&k2) - // Add left-to-right using the NAF optimization. See algorithm 3.77 from - // [GECC]. + // Add left-to-right using the endomorphism and wNAF optimizations. See + // algorithms 3.36 and 3.77 from [GECC]. // // Point Q = ∞ (point at infinity). + // + // Like ordinary binary scalar multiplication, the accumulator is doubled + // once per processed bit position. However, instead of each bit + // contributing either 0 or 1 copies of the point, each non-zero wNAF digit + // contributes a precomputed odd multiple while the preceding doublings + // supply the required power-of-two scaling. + // + // ±P, ±3P, ±5P, ... + // + // Since non-zero digits are sparse and separated by several zero digits, + // relatively few point additions are required compared to the ordinary + // binary representation. var q JacobianPoint - m := k1Len - if m < k2Len { - m = k2Len + maxBits := k1NAF.bits + if maxBits < k2NAF.bits { + maxBits = k2NAF.bits } - for i := 0; i < m; i++ { - // Since k1 and k2 are potentially different lengths and the calculation - // is being done left to right, pad the front of the shorter one with - // 0s. - var k1BytePos, k1ByteNeg, k2BytePos, k2ByteNeg byte - if i >= m-k1Len { - k1BytePos, k1ByteNeg = k1PosNAF[i-(m-k1Len)], k1NegNAF[i-(m-k1Len)] - } - if i >= m-k2Len { - k2BytePos, k2ByteNeg = k2PosNAF[i-(m-k2Len)], k2NegNAF[i-(m-k2Len)] + for bit := maxBits; bit > 0; bit-- { + // Q = 2 * Q + DoubleNonConst(&q, &q) + + // The encoded wNAF digit of k1 serves directly as an index into the + // precomputed odd multiples of P. + if code := k1NAF.codes[bit-1]; code != 0 { + AddNonConst(&q, &pPrecomps[code], &q) } - for mask := uint8(1 << 7); mask > 0; mask >>= 1 { - // Q = 2 * Q - DoubleNonConst(&q, &q) - - // Add or subtract the first point based on the signed digit of the - // NAF representation of k1 at this bit position. - // - // +1: Q = Q + p1 - // -1: Q = Q - p1 - // 0: Q = Q (no change) - if k1BytePos&mask == mask { - AddNonConst(&q, p1, &q) - } else if k1ByteNeg&mask == mask { - AddNonConst(&q, p1Neg, &q) - } - - // Add or subtract the second point based on the signed digit of the - // NAF representation of k2 at this bit position. - // - // +1: Q = Q + p2 - // -1: Q = Q - p2 - // 0: Q = Q (no change) - if k2BytePos&mask == mask { - AddNonConst(&q, p2, &q) - } else if k2ByteNeg&mask == mask { - AddNonConst(&q, p2Neg, &q) - } + // Likewise, the encoded wNAF digit of k2 serves directly as an index + // into the precomputed odd multiples of φ(P). + if code := k2NAF.codes[bit-1]; code != 0 { + AddNonConst(&q, &phiPrecomps[code], &q) } } From 721ed07ea2d081a6f93d093ef726704e0ab632e4 Mon Sep 17 00:00:00 2001 From: Dave Collins Date: Sun, 19 Jul 2026 08:49:14 -0500 Subject: [PATCH 16/16] secp256k1: Remove unused NAF code. This removes the old ordinary NAF implementation that is no longer used now that the implementation has been updated to use a new width-5 wNAF recoder instead. It also removes the associated tests and benchmark accordingly. --- dcrec/secp256k1/curve.go | 136 ---------------------------- dcrec/secp256k1/curve_bench_test.go | 13 --- dcrec/secp256k1/curve_test.go | 119 ------------------------ 3 files changed, 268 deletions(-) diff --git a/dcrec/secp256k1/curve.go b/dcrec/secp256k1/curve.go index 9c2830264..e4cb3f514 100644 --- a/dcrec/secp256k1/curve.go +++ b/dcrec/secp256k1/curve.go @@ -1171,142 +1171,6 @@ func wnaf(k *ModNScalar) wnafScalar { return result } -// nafScalar represents a positive integer up to a maximum value of 2^256 - 1 -// encoded in non-adjacent form. -// -// NAF is a signed-digit representation where each digit can be +1, 0, or -1. -// -// In order to efficiently encode that information, this type uses two arrays, a -// "positive" array where set bits represent the +1 signed digits and a -// "negative" array where set bits represent the -1 signed digits. 0 is -// represented by neither array having a bit set in that position. -// -// The Pos and Neg methods return the aforementioned positive and negative -// arrays, respectively. -type nafScalar struct { - // pos houses the positive portion of the representation. An additional - // byte is required for the positive portion because the NAF encoding can be - // up to 1 bit longer than the normal binary encoding of the value. - // - // neg houses the negative portion of the representation. Even though the - // additional byte is not required for the negative portion, since it can - // never exceed the length of the normal binary encoding of the value, - // keeping the same length for positive and negative portions simplifies - // working with the representation and allows extra conditional branches to - // be avoided. - // - // start and end specify the starting and ending index to use within the pos - // and neg arrays, respectively. This allows fixed size arrays to be used - // versus needing to dynamically allocate space on the heap. - // - // NOTE: The fields are defined in the order that they are to minimize the - // padding on 32-bit and 64-bit platforms. - pos [33]byte - start, end uint8 - neg [33]byte -} - -// Pos returns the bytes of the encoded value with bits set in the positions -// that represent a signed digit of +1. -func (s *nafScalar) Pos() []byte { - return s.pos[s.start:s.end] -} - -// Neg returns the bytes of the encoded value with bits set in the positions -// that represent a signed digit of -1. -func (s *nafScalar) Neg() []byte { - return s.neg[s.start:s.end] -} - -// naf takes a positive integer up to a maximum value of 2^256 - 1 and returns -// its non-adjacent form (NAF), which is a unique signed-digit representation -// such that no two consecutive digits are nonzero. See the documentation for -// the returned type for details on how the representation is encoded -// efficiently and how to interpret it -// -// NAF is useful in that it has the fewest nonzero digits of any signed digit -// representation, only 1/3rd of its digits are nonzero on average, and at least -// half of the digits will be 0. -// -// The aforementioned properties are particularly beneficial for optimizing -// elliptic curve point multiplication because they effectively minimize the -// number of required point additions in exchange for needing to perform a mix -// of fewer point additions and subtractions and possibly one additional point -// doubling. This is an excellent tradeoff because subtraction of points has -// the same computational complexity as addition of points and point doubling is -// faster than both. -func naf(k []byte) nafScalar { - // Strip leading zero bytes. - for len(k) > 0 && k[0] == 0x00 { - k = k[1:] - } - - // The non-adjacent form (NAF) of a positive integer k is an expression - // k = ∑_(i=0, l-1) k_i * 2^i where k_i ∈ {0,±1}, k_(l-1) != 0, and no two - // consecutive digits k_i are nonzero. - // - // The traditional method of computing the NAF of a positive integer is - // given by algorithm 3.30 in [GECC]. It consists of repeatedly dividing k - // by 2 and choosing the remainder so that the quotient (k−r)/2 is even - // which ensures the next NAF digit is 0. This requires log_2(k) steps. - // - // However, in [BRID], Prodinger notes that a closed form expression for the - // NAF representation is the bitwise difference 3k/2 - k/2. This is more - // efficient as it can be computed in O(1) versus the O(log(n)) of the - // traditional approach. - // - // The following code makes use of that formula to compute the NAF more - // efficiently. - // - // To understand the logic here, observe that the only way the NAF has a - // nonzero digit at a given bit is when either 3k/2 or k/2 has a bit set in - // that position, but not both. In other words, the result of a bitwise - // xor. This can be seen simply by considering that when the bits are the - // same, the subtraction is either 0-0 or 1-1, both of which are 0. - // - // Further, observe that the "+1" digits in the result are contributed by - // 3k/2 while the "-1" digits are from k/2. So, they can be determined by - // taking the bitwise and of each respective value with the result of the - // xor which identifies which bits are nonzero. - // - // Using that information, this loops backwards from the least significant - // byte to the most significant byte while performing the aforementioned - // calculations by propagating the potential carry and high order bit from - // the next word during the right shift. - kLen := len(k) - var result nafScalar - var carry uint8 - for byteNum := kLen - 1; byteNum >= 0; byteNum-- { - // Calculate k/2. Notice the carry from the previous word is added and - // the low order bit from the next word is shifted in accordingly. - kc := uint16(k[byteNum]) + uint16(carry) - var nextWord uint8 - if byteNum > 0 { - nextWord = k[byteNum-1] - } - halfK := kc>>1 | uint16(nextWord<<7) - - // Calculate 3k/2 and determine the non-zero digits in the result. - threeHalfK := kc + halfK - nonZeroResultDigits := threeHalfK ^ halfK - - // Determine the signed digits {0, ±1}. - result.pos[byteNum+1] = uint8(threeHalfK & nonZeroResultDigits) - result.neg[byteNum+1] = uint8(halfK & nonZeroResultDigits) - - // Propagate the potential carry from the 3k/2 calculation. - carry = uint8(threeHalfK >> 8) - } - result.pos[0] = carry - - // Set the starting and ending positions within the fixed size arrays to - // identify the bytes that are actually used. This is important since the - // encoding is big endian and thus trailing zero bytes changes its value. - result.start = 1 - carry - result.end = uint8(kLen + 1) - return result -} - // wNAFPrecompTableSize is the size of the wNAF precomputed point table. const wNAFPrecompTableSize = 1 << (wNAFWidth - 1) diff --git a/dcrec/secp256k1/curve_bench_test.go b/dcrec/secp256k1/curve_bench_test.go index 5bcd50e17..7d8bf3089 100644 --- a/dcrec/secp256k1/curve_bench_test.go +++ b/dcrec/secp256k1/curve_bench_test.go @@ -216,19 +216,6 @@ func BenchmarkScalarMultNonConst(b *testing.B) { } } -// BenchmarkNAF benchmarks conversion of a positive integer into its -// non-adjacent form representation. -func BenchmarkNAF(b *testing.B) { - k := fromHex("d74bf844b0862475103d96a611cf2d898447e288d34b360bc885cb8ce7c00575") - kBytes := k.Bytes() - - b.ReportAllocs() - b.ResetTimer() - for i := 0; i < b.N; i++ { - naf(kBytes) - } -} - // BenchmarkWNAF benchmarks conversion of a non-negative integer into its // width-w windowed non-adjacent form representation, where w is [wNAFWidth]. func BenchmarkWNAF(b *testing.B) { diff --git a/dcrec/secp256k1/curve_test.go b/dcrec/secp256k1/curve_test.go index 607e1ce32..c852b0871 100644 --- a/dcrec/secp256k1/curve_test.go +++ b/dcrec/secp256k1/curve_test.go @@ -634,125 +634,6 @@ func TestDoubleJacobian(t *testing.T) { } } -// checkNAFEncoding returns an error if the provided positive and negative -// portions of an overall NAF encoding do not adhere to the requirements or they -// do not sum back to the provided original value. -func checkNAFEncoding(pos, neg []byte, origValue *big.Int) error { - // NAF must not have a leading zero byte and the number of negative - // bytes must not exceed the positive portion. - if len(pos) > 0 && pos[0] == 0 { - return fmt.Errorf("positive has leading zero -- got %x", pos) - } - if len(neg) > len(pos) { - return fmt.Errorf("negative has len %d > pos len %d", len(neg), - len(pos)) - } - - // Ensure the result doesn't have any adjacent non-zero digits. - gotPos := new(big.Int).SetBytes(pos) - gotNeg := new(big.Int).SetBytes(neg) - posOrNeg := new(big.Int).Or(gotPos, gotNeg) - prevBit := posOrNeg.Bit(0) - for bit := 1; bit < posOrNeg.BitLen(); bit++ { - thisBit := posOrNeg.Bit(bit) - if prevBit == 1 && thisBit == 1 { - return fmt.Errorf("adjacent non-zero digits found at bit pos %d", - bit-1) - } - prevBit = thisBit - } - - // Ensure the resulting positive and negative portions of the overall - // NAF representation sum back to the original value. - gotValue := new(big.Int).Sub(gotPos, gotNeg) - if origValue.Cmp(gotValue) != 0 { - return fmt.Errorf("pos-neg is not original value: got %x, want %x", - gotValue, origValue) - } - - return nil -} - -// TestNAF ensures encoding various edge cases and values to non-adjacent form -// produces valid results. -func TestNAF(t *testing.T) { - tests := []struct { - name string // test description - in string // hex encoded test value - }{{ - name: "empty is zero", - in: "", - }, { - name: "zero", - in: "00", - }, { - name: "just before first carry", - in: "aa", - }, { - name: "first carry", - in: "ab", - }, { - name: "leading zeroes", - in: "002f20569b90697ad471c1be6107814f53f47446be298a3a2a6b686b97d35cf9", - }, { - name: "257 bits when NAF encoded", - in: "c000000000000000000000000000000000000000000000000000000000000001", - }, { - name: "32-byte scalar", - in: "6df2b5d30854069ccdec40ae022f5c948936324a4e9ebed8eb82cfd5a6b6d766", - }, { - name: "first term of balanced length-two representation #1", - in: "b776e53fb55f6b006a270d42d64ec2b1", - }, { - name: "second term balanced length-two representation #1", - in: "d6cc32c857f1174b604eefc544f0c7f7", - }, { - name: "first term of balanced length-two representation #2", - in: "45c53aa1bb56fcd68c011e2dad6758e4", - }, { - name: "second term of balanced length-two representation #2", - in: "a2e79d200f27f2360fba57619936159b", - }} - - for _, test := range tests { - // Ensure the resulting positive and negative portions of the overall - // NAF representation adhere to the requirements of NAF encoding and - // they sum back to the original value. - result := naf(hexToBytes(test.in)) - pos, neg := result.Pos(), result.Neg() - if err := checkNAFEncoding(pos, neg, fromHex(test.in)); err != nil { - t.Errorf("%q: %v", test.name, err) - } - } -} - -// TestNAFRandom ensures that encoding randomly-generated values to non-adjacent -// form produces valid results. -func TestNAFRandom(t *testing.T) { - // Use a unique random seed each test instance and log it if the tests fail. - seed := time.Now().Unix() - rng := mrand.New(mrand.NewSource(seed)) - defer func(t *testing.T, seed int64) { - if t.Failed() { - t.Logf("random seed: %d", seed) - } - }(t, seed) - - for i := 0; i < 100; i++ { - // Ensure the resulting positive and negative portions of the overall - // NAF representation adhere to the requirements of NAF encoding and - // they sum back to the original value. - bigIntVal, modNVal := randIntAndModNScalar(t, rng) - valBytes := modNVal.Bytes() - result := naf(valBytes[:]) - pos, neg := result.Pos(), result.Neg() - if err := checkNAFEncoding(pos, neg, bigIntVal); err != nil { - t.Fatalf("encoding err: %v\nin: %x\npos: %x\nneg: %x", err, - bigIntVal, pos, neg) - } - } -} - // decodeWNAFDigit converts the provided wNAF code to the signed digit it // represents. func decodeWNAFDigit(code uint8) int8 {