fieldmul.gno
16.68 Kb · 378 lines
1package bls254
2
3// Field multiplication mod P. TWO implementations, and this header records which
4// one the drand path uses and the measurement that decided it. PHASE 2 SETTLED
5// THIS; please do not re-litigate it from a gas table that came from somewhere
6// else, because that is exactly how it went wrong the first time.
7//
8// The rest of this package gets away with crypto/modexp for every field
9// operation it needs, because try-and-increment hash-to-G1 only ever needs
10// powers (x^3, sqrt) and one addition. The RFC 9380 SVDW map that the drand
11// path needs (see drandcurve.gno) is different: it needs about a dozen general
12// products a*b per field element mapped, and modexp cannot multiply two
13// distinct operands.
14//
15// Two ways out:
16//
17// 1. The identity a*b == ((a+b)^2 - (a-b)^2) / 4, which crypto/modexp can do
18// because it can square. Two native modexp calls, one modular subtraction and
19// a division by 4, which mod P is two halvings and not a third multiply.
20// Implemented below as mulModPIdentity / mulModPLimbs.
21// 2. A hand-rolled 4-limb Montgomery multiply over math/bits.Mul64. No native
22// call, 32 64x64 multiplies per montMul. Implemented below as mulModP.
23//
24// (1) WINS, AND IT IS WHAT drandcurve.gno's fpMul AND fpSqr BIND TO.
25//
26// Measured on this toolchain (/home/johnny/go/bin/gno, gno v1.2.0 stdlibs) from
27// gno test's own per-test GAS line, 100 iterations per test with the 349,766-gas
28// empty-test baseline subtracted. The tests are in drandgas_test.gno; re-run them
29// rather than trusting these numbers.
30//
31// operation gas each
32// mulModP 4-limb Montgomery 2,200,369
33// mulModPIdentity identity, byte domain 519,791 <- 4.23x cheaper
34// mulModPLimbs identity, limb domain 387,943
35// sqrModP Montgomery a*a 2,200,184
36// sqrModPNative one native modexp 69,602 <- 31.6x cheaper
37// modexp.ModExp(x, 2, P) bare 23,145
38// modexp.ModExp(x, (P+1)/4, P) bare 134,372
39// toLimbs / fromLimbs 39,328 / 42,308
40// halveModP 21,035
41//
42// End to end, net of the same baseline:
43//
44// Montgomery identity factor
45// VerifyDrandOnce 63,846,571 13,396,451 4.77x
46// VerifyDrandFour / 4 68,401,101 15,106,773 4.53x
47// HashToG1DrandOnce 62,394,145 12,004,207 5.20x
48// MapToPointSVDWOnce 35,797,798 6,657,484 5.38x
49// r/splash/beacon Publish 74,891,351 16,671,000* 4.49x
50//
51// THE IDENTITY COLUMN ROSE ON 2026-09-29 AND THE REASON IS ONE LINE OF CODE, NOT DRIFT.
52// mulModPLimbs gained the input guard AUDIT-2026-09-28.md finding 21 asked for -- two
53// cmp calls against pLimbs, measured at 8,500 gas the pair -- so every multiply costs
54// 8,500 more and the end-to-end figures rise by 8,500 times the number of multiplies
55// they perform (16 in VerifyDrandOnce, 9 in the SVDW map). VerifyDrand is up 1.0 %.
56// The previous figures were 511,291 / 379,443 / 13,260,451 / 14,964,398 / 11,868,207 /
57// 6,580,984, and they are recorded here rather than silently replaced so that a reader
58// who remembers them learns why they moved. mulModP, sqrModP and sqrModPNative are
59// unchanged because the guard is not on their path.
60//
61// *THE Publish ROW IS DERIVED, NOT RE-ISOLATED, and it is the only number in this file
62// that is. A Publish is its own fixed cost plus exactly one VerifyDrand, so it rises by
63// the same 136,000; r/splash/beacon's own suite has no single-Publish gas test to read
64// it off, and this change may not touch that realm. Re-isolate it there before quoting
65// it as measured.
66//
67// WHY THE ORIGINAL CHOICE WENT THE OTHER WAY. (2) was specified on the strength of
68// a gas table reading "modexp: Base 58000 + 24647680*len(mod)/1024", which would
69// have made a modexp call cost ~800k. That table is from a different gno revision
70// -- this host has two pinned toolchains, and it was read off GnoBuilder's
71// c4c72fdd instead of SPLASH's v1.2.0. On v1.2.0 the crypto/modexp row is Base
72// 1400, Slope 2200 on SizeModExpWork (a count of modular multiplications) and
73// Slope2 3500 per 1024 bytes of modulus: a 254-bit squaring costs 23k, not 800k.
74// Interpreted 64-bit limb arithmetic is the expensive thing here, not the native
75// call. If you are about to change this back, measure on the toolchain you will
76// deploy with first.
77//
78// MONTGOMERY IS KEPT, NOT DEAD CODE. mulModP is the independent oracle that
79// TestMulModPIdentityMatchesMontgomery compares the identity multiply against:
80// two implementations sharing no arithmetic, checked byte-for-byte over every
81// ordered pair of 15 limb-boundary edge values, the python-derived vector table,
82// 120 in-VM pseudorandom pairs and all 11 rejection cases, with the tally
83// asserting that a+b wrapping, a-b going negative and both halvings' odd and even
84// branches were each actually reached. Deleting mulModP would delete that oracle.
85//
86// WHAT IS STILL ON THE TABLE, measured rather than guessed. The multiply is now
87// dominated by byte<->limb marshalling, not arithmetic: mulModPLimbs is 387,943
88// against mulModPIdentity's 519,791, purely because the caller holds limbs. Moving
89// mapToPointSVDW itself into the limb domain would save ~132k on each of its ~8
90// products and ~120k on each of its ~6 additions, but would cost ~34k on each of
91// its ~5 squarings (a squaring is cheapest in the byte domain, where it is one
92// bare modexp) and ~42-81k on each of its 3 full exponentiations. Net estimate
93// ~1.4M off a 6.66M map, so VerifyDrand ~13.40M -> ~10.5M: a further 21% for a
94// rewrite of the most correctness-critical function in the package. Not taken.
95//
96// AND SUCH A REWRITE WOULD BE THE FIRST REAL LIMB-DOMAIN CALLER OF mulModPLimbs, which
97// is precisely who finding 21's guard is for: the byte-domain wrappers would no longer
98// be standing between the SVDW map and the unreduced-operand case.
99// 0.0165 GNOT a relay is already 4.5x inside the budget that prompted this work.
100//
101// One cheap, safe win was left on the table deliberately: inlining mac into the
102// two loops below measured 360k instead of 503k gas for mul512, a 28% saving.
103// mul512 is now off the hot path entirely -- it is only reached through the
104// oracle -- so readability of the reduction loop wins outright and mac stays a
105// function.
106//
107// Everything here is unexported. mulModP, mulModPIdentity and sqrModPNative all
108// take and return 32-byte big-endian slices, matching modexp and the point
109// encodings used everywhere else in this package, and all three return nil for an
110// operand that is not exactly 32 bytes or is >= P.
111
112import (
113 "crypto/modexp"
114 "math/bits"
115)
116
117const (
118 // pInvNeg64 is -P^{-1} mod 2^64, the Montgomery reduction constant.
119 //
120 // Derived, not recalled:
121 // python3 -c "p=0x30644e72e131a029b85045b68181585d97816a916871ca8d3c208c16d87cfd47;
122 // print(hex((-pow(p,-1,1<<256)) % (1<<64)))"
123 // -> 0x87d20782e4866389
124 // TestMontgomeryConstants re-derives it inside the VM from pLimbs alone,
125 // so a typo here fails the test suite rather than silently producing
126 // wrong products.
127 pInvNeg64 = 0x87d20782e4866389
128
129 // r2Hex is R^2 mod P where R = 2^256, i.e. 2^512 mod P. montMul(x, r2)
130 // converts x into Montgomery form; montMul(x, 1) converts back.
131 //
132 // python3 -c "p=0x30644e72e131a029b85045b68181585d97816a916871ca8d3c208c16d87cfd47;
133 // print('%064x' % pow(2,512,p))"
134 r2Hex = "06d89f71cab8351f47ab1eff0a417ff6b5e71911d44501fbf32cfc5b538afa89"
135)
136
137var r2Limbs = toLimbs(mustDecodeHex(r2Hex))
138
139// mac returns the 128-bit a*b + c + carryIn as (hi, lo).
140//
141// It cannot overflow: a*b <= (2^64-1)^2 = 2^128 - 2^65 + 1, so hi <= 2^64 - 2
142// and the two carries that can come out of the low word still leave hi below
143// 2^64.
144func mac(a, b, c, carryIn uint64) (hi, lo uint64) {
145 hi, lo = bits.Mul64(a, b)
146 var c1 uint64
147 lo, c1 = bits.Add64(lo, c, 0)
148 hi, _ = bits.Add64(hi, 0, c1)
149 lo, c1 = bits.Add64(lo, carryIn, 0)
150 hi, _ = bits.Add64(hi, 0, c1)
151 return hi, lo
152}
153
154// mul512 is the schoolbook 4x4 -> 8 limb product. Limbs are little-endian
155// (index 0 least significant), matching toLimbs.
156func mul512(a, b [4]uint64) [8]uint64 {
157 var t [8]uint64
158 for i := 0; i < 4; i++ {
159 var carry uint64
160 for j := 0; j < 4; j++ {
161 carry, t[i+j] = mac(a[j], b[i], t[i+j], carry)
162 }
163 // t[i+4] has not been written by any earlier i, so this is an
164 // assignment rather than an add.
165 t[i+4] = carry
166 }
167 return t
168}
169
170// montReduce returns t * R^{-1} mod P for t < P*R, by separated-operand-scanning
171// Montgomery reduction.
172func montReduce(t [8]uint64) [4]uint64 {
173 // hi9 is the notional ninth limb. The loop can carry past t[7] because
174 // t + m_i*P*2^(64i) can exceed 2^512 even though the final quotient does
175 // not; see the note on the return value below.
176 var hi9 uint64
177
178 for i := 0; i < 4; i++ {
179 m := t[i] * pInvNeg64
180 var carry uint64
181 for j := 0; j < 4; j++ {
182 carry, t[i+j] = mac(m, pLimbs[j], t[i+j], carry)
183 }
184 // t[i] is now zero by construction. Propagate the high carry.
185 for k := i + 4; k < 8 && carry != 0; k++ {
186 t[k], carry = bits.Add64(t[k], carry, 0)
187 }
188 hi9 += carry
189 }
190
191 r := [4]uint64{t[4], t[5], t[6], t[7]}
192
193 // For a, b < P the input to this function is a*b < P^2 < P*R, so the
194 // quotient is below P + P^2/R < 2P. P < 2^254, so 2P < 2^255 and hi9 is
195 // necessarily zero. Fold it into the reduction condition anyway: if a
196 // caller ever feeds an unreduced operand the answer stays in range for
197 // the single subtraction below instead of silently losing 2^256.
198 if hi9 != 0 || cmp(r, pLimbs) >= 0 {
199 r = sub(r, pLimbs)
200 }
201 return r
202}
203
204// montMul returns a*b*R^{-1} mod P for a, b reduced mod P.
205func montMul(a, b [4]uint64) [4]uint64 {
206 return montReduce(mul512(a, b))
207}
208
209// mulModP returns (a*b) mod P. a and b must be 32-byte big-endian values below
210// P; the result is 32 bytes big-endian. Returns nil if either input is not
211// exactly 32 bytes or is >= P -- an unreduced operand is a caller bug, and
212// returning a plausible-looking wrong product for it is the one outcome this
213// package cannot afford.
214func mulModP(a, b []byte) []byte {
215 if len(a) != 32 || len(b) != 32 {
216 return nil
217 }
218 al := toLimbs(a)
219 bl := toLimbs(b)
220 if cmp(al, pLimbs) >= 0 || cmp(bl, pLimbs) >= 0 {
221 return nil
222 }
223 // montMul(a, b) = a*b*R^{-1}; multiplying that by R^2 gives a*b.
224 return fromLimbs(montMul(montMul(al, bl), r2Limbs))
225}
226
227// sqrModP returns (a*a) mod P, with the same contract as mulModP.
228func sqrModP(a []byte) []byte { return mulModP(a, a) }
229
230// ---- Approach (1): the squaring identity, via native crypto/modexp ----
231//
232// This is the alternative the header calls the easier one. It is here because it
233// is also, measured on this toolchain, the much cheaper one -- see the table at
234// the bottom of this file.
235//
236// a*b == ((a+b)^2 - (a-b)^2) / 4 (mod P)
237//
238// The identity is exact over the integers: (a+b)^2 - (a-b)^2 = 4ab. Reducing
239// a+b and a-b mod P first does not change their squares mod P, so the numerator
240// is congruent to 4ab and the division by 4 -- which mod P means multiplying by
241// the inverse of 4, i.e. halving twice -- recovers ab.
242
243// halveModP returns (x / 2) mod P, i.e. x * inverse(2) mod P, for x reduced mod P.
244//
245// P is odd, so inverse(2) exists. If x is even the integer x/2 is already the
246// answer. If x is odd then x+P is even and (x+P)/2 is the answer, because
247// 2*((x+P)/2) = x+P == x (mod P).
248//
249// No overflow: P < 2^254 and x < P, so x+P < 2^255. The carry out of add is
250// therefore always zero; it is folded into the shift anyway so that an unreduced
251// input cannot silently lose the top bit.
252func halveModP(x [4]uint64) [4]uint64 {
253 var carry uint64
254 if x[0]&1 == 1 {
255 x, carry = add(x, pLimbs)
256 }
257 // Shift right by one, pulling in carry as bit 256.
258 var r [4]uint64
259 r[0] = x[0]>>1 | x[1]<<63
260 r[1] = x[1]>>1 | x[2]<<63
261 r[2] = x[2]>>1 | x[3]<<63
262 r[3] = x[3]>>1 | carry<<63
263 return r
264}
265
266// sqrModPNative returns (a*a) mod P using one native modexp, with exactly
267// mulModP's input contract: nil unless a is 32 bytes and below P.
268//
269// modexp.ModExp returns len(modulus) == 32 bytes, left-padded, so the result
270// needs no adjustment.
271func sqrModPNative(a []byte) []byte {
272 if len(a) != 32 {
273 return nil
274 }
275 if cmp(toLimbs(a), pLimbs) >= 0 {
276 return nil
277 }
278 return modexp.ModExp(a, twoBytes, pBytes)
279}
280
281// twoBytes is the exponent 2. A one-byte exponent is what modexp's gas row
282// charges least for.
283var twoBytes = []byte{2}
284
285// mulModPIdentity returns (a*b) mod P via the squaring identity above. Its
286// contract is identical to mulModP's, including returning nil for a non-32-byte
287// or unreduced operand: TestMulModPIdentityMatchesMontgomery pins the two
288// functions to agreeing on every input, rejections included.
289func mulModPIdentity(a, b []byte) []byte {
290 if len(a) != 32 || len(b) != 32 {
291 return nil
292 }
293 al := toLimbs(a)
294 bl := toLimbs(b)
295 if cmp(al, pLimbs) >= 0 || cmp(bl, pLimbs) >= 0 {
296 return nil
297 }
298
299 return fromLimbs(mulModPLimbs(al, bl))
300}
301
302// mulModPLimbs is the identity multiply in the limb domain, for callers that
303// already hold reduced limbs and want reduced limbs back. Both operands MUST be
304// below P, AND SINCE 2026-09-29 THAT IS CHECKED HERE RATHER THAN ASSUMED.
305//
306// THE DEFECT THAT PUT THE CHECK IN, AUDIT-2026-09-28.md finding 21. This comment used
307// to read "there is no check, because the byte-domain wrappers do it and the SVDW map
308// only ever feeds it values it has already reduced". Both halves of that were true, and
309// the function was still THE ONLY PRODUCER IN THIS PACKAGE THAT WAS NOT FAIL-CLOSED.
310// Fed an unreduced b with b > a + P, subModP's `a < b` branch computes a + (P - b), and
311// P - b wraps mod 2^256 because b > P -- so the result is NOT congruent to a - b, the
312// numerator is not congruent to 4ab, and the function returns A PLAUSIBLE-LOOKING WRONG
313// PRODUCT. This package's own header calls that the one outcome it cannot afford: every
314// other entry point returns nil rather than a value that is merely probably right.
315// mulModP, mulModPIdentity and sqrModPNative all check; this did not.
316//
317// IT IS A PANIC AND NOT A nil BECAUSE THE RETURN TYPE HAS NO nil. [4]uint64 has no
318// invalid value to hand back, and its zero value is a legitimate product, so the only
319// fail-closed answer available is to stop. That is the choice hashToFieldDrand and
320// hash_to_g1 already make for their own impossible states.
321//
322// IT IS UNREACHABLE TODAY AND THE GUARD IS THEREFORE A CONTINGENCY, NOT A FIX TO A LIVE
323// BUG: the only non-test caller is mulModPIdentity, which performs exactly this
324// comparison on the same values first, and drandcurve.gno's fpMul and fpSqr still bind
325// to the input-checked mulModPIdentity and sqrModPNative. NOTHING ON THE LIVE PATH
326// CHANGED, and the measured cost of the redundant pair of cmp calls is in the header's
327// table. What changed is that a FUTURE caller in the limb domain -- which is the whole
328// reason this function is factored out at all -- cannot silently get a wrong answer.
329//
330// This exists because the byte<->limb conversions, not the arithmetic, are what
331// this multiply actually costs. toLimbs and fromLimbs measure ~53k and ~58k gas
332// each (see drandgas_test.gno); the byte-domain wrapper pays five conversions,
333// this core pays the four that the modexp boundary makes unavoidable.
334func mulModPLimbs(al, bl [4]uint64) [4]uint64 {
335 if cmp(al, pLimbs) >= 0 || cmp(bl, pLimbs) >= 0 {
336 panic("bls254: mulModPLimbs: operand is not reduced mod P")
337 }
338
339 // s = (a+b) mod P, d = (a-b) mod P. addModP handles the carry out of a+b
340 // exceeding P; subModP handles a < b by adding P back.
341 s := fromLimbs(addModP(al, bl, pLimbs))
342 d := fromLimbs(subModP(al, bl))
343
344 // Two native squarings. Both operands are already reduced.
345 s2 := modexp.ModExp(s, twoBytes, pBytes)
346 d2 := modexp.ModExp(d, twoBytes, pBytes)
347
348 // n == 4ab (mod P); halve twice to divide by 4.
349 n := subModP(toLimbs(s2), toLimbs(d2))
350 return halveModP(halveModP(n))
351}
352
353// subModP returns (a - b) mod P for a, b reduced mod P.
354//
355// THE PRECONDITION IS LOAD-BEARING AND IS WHERE FINDING 21 ACTUALLY LIVED. On the
356// `a < b` branch it computes a + (P - b), and sub() wraps mod 2^256: so for b > P the
357// term P - b is a huge wrong number and the result is not congruent to a - b. There is
358// no way to detect that from inside this function -- every [4]uint64 is a valid return
359// value -- so the contract is enforced by its callers, and mulModPLimbs states it as a
360// panic rather than assuming it.
361func subModP(a, b [4]uint64) [4]uint64 {
362 if cmp(a, b) >= 0 {
363 return sub(a, b)
364 }
365 // a < b: a - b + P, computed as a + (P - b). P - b is in [1, P] and a is
366 // below b <= P-1, so the sum is below P and no carry escapes.
367 d, _ := add(a, sub(pLimbs, b))
368 return d
369}
370
371// negModP returns (-a) mod P for a reduced mod P.
372func negModP(a [4]uint64) [4]uint64 {
373 var zero [4]uint64
374 if cmp(a, zero) == 0 {
375 return zero
376 }
377 return sub(pLimbs, a)
378}
