nautilus_lighter/signing/field/
goldilocks.rs

1// -------------------------------------------------------------------------------------------------
2//  Copyright (C) 2015-2026 Nautech Systems Pty Ltd. All rights reserved.
3//  https://nautechsystems.io
4//
5//  Licensed under the GNU Lesser General Public License Version 3.0 (the "License");
6//  You may not use this file except in compliance with the License.
7//  You may obtain a copy of the License at https://www.gnu.org/licenses/lgpl-3.0.en.html
8//
9//  Unless required by applicable law or agreed to in writing, software
10//  distributed under the License is distributed on an "AS IS" BASIS,
11//  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12//  See the License for the specific language governing permissions and
13//  limitations under the License.
14// -------------------------------------------------------------------------------------------------
15
16//! Goldilocks prime field `Fp = GF(p)`, with `p = 2^64 - 2^32 + 1`.
17//!
18//! Elements are stored in Montgomery form internally so multiplication reduces
19//! to a single 64x64 -> 128-bit multiply followed by a fixed-shape Montgomery
20//! reduction; values in the canonical `0..p-1` range are never observed before
21//! a deliberate `to_u64`/`to_le_bytes` call. The arithmetic core (`+`, `-`,
22//! `*`, `neg`, `square`, `msquare`, `pow`, `invert`) executes as a
23//! straight-line sequence of arithmetic and bitwise ops with no data-dependent
24//! branches, so timing leaks no information about field operands. The
25//! Tonelli-Shanks [`Fp::sqrt`] is variable-time over its input and is only
26//! consumed by [`super::Fp5::sqrt`] / [`super::Fp5::canonical_sqrt`], which
27//! the curve `Point::decode` calls on public-input `w` values; secret-input
28//! sqrt is not part of the signing critical path.
29
30use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
31
32/// Goldilocks prime modulus: `p = 2^64 - 2^32 + 1`.
33pub const MODULUS: u64 = 0xFFFF_FFFF_0000_0001;
34
35/// `R^2 mod p` with `R = 2^64`. Used to lift a `u64` into Montgomery form.
36const R2: u64 = 0xFFFF_FFFE_0000_0001;
37
38/// 2-adicity of `p - 1`: `p - 1 = 2^32 * (2^32 - 1)`.
39const TWO_ADICITY: u32 = 32;
40
41/// Generator of the unique subgroup of order `2^32` in `Fp^*`. Used as the
42/// Tonelli-Shanks "non-residue" `z`. Matches the Plonky2 / Lighter convention.
43const POWER_OF_TWO_GENERATOR: u64 = 7_277_203_076_849_721_926;
44
45/// An element of the Goldilocks field `Fp = GF(p)`.
46///
47/// The wrapped `u64` holds the value in Montgomery form (`x * 2^64 mod p`),
48/// always reduced into `0..p-1`. Two `Fp` instances are equal iff their
49/// Montgomery limbs are equal, so `PartialEq` is correct and constant-time.
50#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
51pub struct Fp(pub(super) u64);
52
53impl Fp {
54    /// Additive identity.
55    pub const ZERO: Self = Self::from_u64_reduce(0);
56
57    /// Multiplicative identity.
58    pub const ONE: Self = Self::from_u64_reduce(1);
59
60    /// Element `-1 mod p`.
61    pub const MINUS_ONE: Self = Self::from_u64_reduce(MODULUS - 1);
62
63    /// Montgomery reduction: given `x` with `x < p * 2^64`, return `x * 2^-64 mod p`
64    /// in canonical `0..p-1` form.
65    #[inline(always)]
66    const fn montyred(x: u128) -> u64 {
67        let xl = x as u64;
68        let xh = (x >> 64) as u64;
69        let (a, e) = xl.overflowing_add(xl << 32);
70        let b = a.wrapping_sub(a >> 32).wrapping_sub(e as u64);
71        let (r, c) = xh.overflowing_sub(b);
72        r.wrapping_sub(0u32.wrapping_sub(c as u32) as u64)
73    }
74
75    #[inline(always)]
76    pub(super) const fn from_montgomery_products(lo: u128, hi: u128) -> Self {
77        // Fold each product's high half with 2^64 = 2^32 - 1 (mod p)
78        Self(Self::montyred(lo + (hi << 32) - hi))
79    }
80
81    /// Build an element from a `u64`, reducing modulo `p`.
82    #[inline(always)]
83    pub const fn from_u64_reduce(v: u64) -> Self {
84        Self(Self::montyred((v as u128) * (R2 as u128)))
85    }
86
87    /// Build an element from an already-canonical `u64` (`v < p`); returns `None` otherwise.
88    #[inline(always)]
89    pub fn from_u64_canonical(v: u64) -> Option<Self> {
90        if v < MODULUS {
91            Some(Self::from_u64_reduce(v))
92        } else {
93            None
94        }
95    }
96
97    /// Return the canonical `u64` representation in `0..p-1`.
98    #[inline(always)]
99    pub const fn to_u64(self) -> u64 {
100        Self::montyred(self.0 as u128)
101    }
102
103    /// Decode an element from 8 little-endian bytes.
104    ///
105    /// Returns `None` if the encoded integer is not in canonical range (`>= p`).
106    #[inline]
107    pub fn try_from_le_bytes(bytes: [u8; 8]) -> Option<Self> {
108        Self::from_u64_canonical(u64::from_le_bytes(bytes))
109    }
110
111    /// Canonical 8-byte little-endian encoding.
112    #[inline]
113    pub fn to_le_bytes(self) -> [u8; 8] {
114        self.to_u64().to_le_bytes()
115    }
116
117    /// Test whether the element is zero.
118    #[inline(always)]
119    pub const fn is_zero(self) -> bool {
120        self.0 == 0
121    }
122
123    /// Constant-time equality: returns `0xFFFF_FFFF_FFFF_FFFF` on equality, `0` otherwise.
124    #[inline(always)]
125    pub const fn ct_eq(self, rhs: Self) -> u64 {
126        let t = self.0 ^ rhs.0;
127        !((((t | t.wrapping_neg()) as i64) >> 63) as u64)
128    }
129
130    /// Branch-free select: returns `a` when `mask == 0` and `b` when
131    /// `mask == u64::MAX`. Behaviour for any other mask value is unspecified;
132    /// the secret-scalar curve primitives only ever pass full-bit masks.
133    #[inline(always)]
134    #[must_use]
135    pub const fn ct_select(mask: u64, a: Self, b: Self) -> Self {
136        Self(a.0 ^ (mask & (a.0 ^ b.0)))
137    }
138
139    #[inline(always)]
140    const fn add_inner(self, rhs: Self) -> Self {
141        let (x1, c1) = self.0.overflowing_sub(MODULUS - rhs.0);
142        let adj = 0u32.wrapping_sub(c1 as u32);
143        Self(x1.wrapping_sub(adj as u64))
144    }
145
146    #[inline(always)]
147    const fn sub_inner(self, rhs: Self) -> Self {
148        let (x1, c1) = self.0.overflowing_sub(rhs.0);
149        let adj = 0u32.wrapping_sub(c1 as u32);
150        Self(x1.wrapping_sub(adj as u64))
151    }
152
153    #[inline(always)]
154    const fn neg_inner(self) -> Self {
155        Self::ZERO.sub_inner(self)
156    }
157
158    #[inline(always)]
159    const fn mul_inner(self, rhs: Self) -> Self {
160        Self(Self::montyred((self.0 as u128) * (rhs.0 as u128)))
161    }
162
163    /// Squaring in `Fp`.
164    #[inline(always)]
165    #[must_use]
166    pub const fn square(self) -> Self {
167        self.mul_inner(self)
168    }
169
170    /// Repeated squaring: returns `self^(2^n)`.
171    #[inline]
172    #[must_use]
173    pub fn msquare(self, n: u32) -> Self {
174        let mut x = self;
175        for _ in 0..n {
176            x = x.square();
177        }
178        x
179    }
180
181    /// Multiplicative inverse via Fermat's little theorem: `x^(p-2)`.
182    ///
183    /// Returns `Fp::ZERO` when called on zero (no panic), matching the
184    /// "inverse-or-zero" convention used by the upstream reference impls.
185    #[must_use]
186    pub fn invert(self) -> Self {
187        // p - 2 = 0xFFFFFFFEFFFFFFFF; addition chain reaches the exponent in 11 mults
188        // and 64 squarings. `xj` denotes `x^(2^j - 1)` at each step.
189        let x = self;
190        let x2 = x * x.square();
191        let x4 = x2 * x2.msquare(2);
192        let x5 = x * x4.square();
193        let x10 = x5 * x5.msquare(5);
194        let x15 = x5 * x10.msquare(5);
195        let x16 = x * x15.square();
196        let x31 = x15 * x16.msquare(15);
197        let x32 = x * x31.square();
198        x32 * x31.msquare(33)
199    }
200
201    /// Exponentiation by an unsigned 64-bit integer, via right-to-left square-and-multiply.
202    #[must_use]
203    pub fn pow(self, mut exp: u64) -> Self {
204        let mut result = Self::ONE;
205        let mut base = self;
206
207        while exp != 0 {
208            if exp & 1 == 1 {
209                result *= base;
210            }
211            base = base.square();
212            exp >>= 1;
213        }
214        result
215    }
216
217    /// Square root in `Fp`.
218    ///
219    /// Returns `Some(s)` such that `s^2 == self` when one exists (with the zero
220    /// element returning `Some(Self::ZERO)`); returns `None` for non-squares.
221    /// Picks one of the two roots: callers wanting a fixed sign must apply
222    /// their own normalization on top.
223    ///
224    /// Implementation is Tonelli-Shanks specialized to the Goldilocks
225    /// `p - 1 = 2^32 * (2^32 - 1)` factorization, with the precomputed
226    /// 2^32-th root-of-unity generator `POWER_OF_TWO_GENERATOR` standing in
227    /// for `z`.
228    #[must_use]
229    pub fn sqrt(self) -> Option<Self> {
230        if self.is_zero() {
231            return Some(Self::ZERO);
232        }
233
234        // Euler's criterion: `self^((p-1)/2)` is `+1` iff `self` is a square.
235        let qr = self.pow((MODULUS - 1) >> 1);
236        if qr == Self::MINUS_ONE {
237            return None;
238        }
239        debug_assert_eq!(qr, Self::ONE);
240
241        let t: u64 = (1u64 << (64 - TWO_ADICITY)) - 1;
242        let mut z = Self::from_u64_reduce(POWER_OF_TWO_GENERATOR);
243        let mut w = self.pow((t - 1) >> 1);
244        let mut x = self * w;
245        let mut b = x * w;
246        let mut v = TWO_ADICITY;
247
248        while b != Self::ONE {
249            let mut k = 0u32;
250            let mut b2k = b;
251
252            while b2k != Self::ONE {
253                b2k = b2k.square();
254                k += 1;
255            }
256
257            let j = v - k - 1;
258            w = z.msquare(j);
259            z = w.square();
260            b *= z;
261            x *= w;
262            v = k;
263        }
264
265        Some(x)
266    }
267}
268
269impl Default for Fp {
270    #[inline]
271    fn default() -> Self {
272        Self::ZERO
273    }
274}
275
276impl Add for Fp {
277    type Output = Self;
278    #[inline(always)]
279    fn add(self, rhs: Self) -> Self {
280        self.add_inner(rhs)
281    }
282}
283
284impl AddAssign for Fp {
285    #[inline(always)]
286    fn add_assign(&mut self, rhs: Self) {
287        *self = self.add_inner(rhs);
288    }
289}
290
291impl Sub for Fp {
292    type Output = Self;
293    #[inline(always)]
294    fn sub(self, rhs: Self) -> Self {
295        self.sub_inner(rhs)
296    }
297}
298
299impl SubAssign for Fp {
300    #[inline(always)]
301    fn sub_assign(&mut self, rhs: Self) {
302        *self = self.sub_inner(rhs);
303    }
304}
305
306impl Neg for Fp {
307    type Output = Self;
308    #[inline(always)]
309    fn neg(self) -> Self {
310        self.neg_inner()
311    }
312}
313
314impl Mul for Fp {
315    type Output = Self;
316    #[inline(always)]
317    fn mul(self, rhs: Self) -> Self {
318        self.mul_inner(rhs)
319    }
320}
321
322impl MulAssign for Fp {
323    #[inline(always)]
324    fn mul_assign(&mut self, rhs: Self) {
325        *self = self.mul_inner(rhs);
326    }
327}
328
329#[cfg(test)]
330mod tests {
331    use proptest::prelude::*;
332    use rstest::rstest;
333    use serde::Deserialize;
334
335    use super::*;
336    use crate::signing::fixtures::{arb_fp, arb_fp_nonzero, hex_to_bytes};
337
338    const VECTORS_JSON: &str = include_str!(concat!(
339        env!("CARGO_MANIFEST_DIR"),
340        "/test_data/signing_field_goldilocks_vectors.json",
341    ));
342
343    #[derive(Debug, Deserialize)]
344    struct Vectors {
345        vectors: Vec<Vector>,
346    }
347
348    #[derive(Debug, Deserialize)]
349    struct Vector {
350        a: String,
351        b: String,
352        e: String,
353        add: String,
354        sub: String,
355        mul: String,
356        neg_a: String,
357        inv_a: String,
358        pow_a_e: String,
359        a_eq_b: bool,
360    }
361
362    fn decode_le8(hex: &str) -> [u8; 8] {
363        let bytes = hex_to_bytes(hex);
364        assert_eq!(bytes.len(), 8, "expected 8 bytes, was {}", bytes.len());
365        let mut out = [0u8; 8];
366        out.copy_from_slice(&bytes);
367        out
368    }
369
370    fn parse_u64(s: &str) -> u64 {
371        if let Some(stripped) = s.strip_prefix("0x") {
372            u64::from_str_radix(stripped, 16).unwrap()
373        } else {
374            s.parse::<u64>().unwrap()
375        }
376    }
377
378    #[rstest]
379    fn modulus_constant_is_goldilocks_prime() {
380        assert_eq!(u128::from(MODULUS), (1u128 << 64) - (1u128 << 32) + 1);
381    }
382
383    #[rstest]
384    fn round_trip_le_bytes_canonical() {
385        for v in [0u64, 1, 42, MODULUS - 1] {
386            let f = Fp::from_u64_canonical(v).unwrap();
387            assert_eq!(f.to_u64(), v);
388            let bytes = f.to_le_bytes();
389            assert_eq!(Fp::try_from_le_bytes(bytes).unwrap(), f);
390        }
391    }
392
393    #[rstest]
394    fn rejects_non_canonical_decoding() {
395        let bad = MODULUS.to_le_bytes();
396        assert!(Fp::try_from_le_bytes(bad).is_none());
397        let worse = u64::MAX.to_le_bytes();
398        assert!(Fp::try_from_le_bytes(worse).is_none());
399    }
400
401    #[rstest]
402    fn invert_zero_returns_zero() {
403        assert_eq!(Fp::ZERO.invert(), Fp::ZERO);
404    }
405
406    #[rstest]
407    fn sqrt_round_trip_for_known_squares() {
408        for v in [1u64, 2, 4, 9, 16, 100, 1_000_000] {
409            let x = Fp::from_u64_reduce(v);
410            let xs = x.square();
411            let s = xs.sqrt().expect("known squares are residues");
412            assert_eq!(s.square(), xs);
413        }
414    }
415
416    #[rstest]
417    fn sqrt_zero_returns_zero() {
418        assert_eq!(Fp::ZERO.sqrt(), Some(Fp::ZERO));
419    }
420
421    #[rstest]
422    fn sqrt_returns_none_for_non_square() {
423        // Construct a guaranteed non-residue by stepping through small candidates
424        // until Euler's criterion fails. Fp has roughly half non-residues, so this
425        // resolves on the first or second probe.
426        let mut v = 2u64;
427
428        loop {
429            let x = Fp::from_u64_reduce(v);
430            if x.pow((MODULUS - 1) >> 1) == Fp::MINUS_ONE {
431                assert_eq!(x.sqrt(), None);
432                break;
433            }
434            v += 1;
435        }
436    }
437
438    #[rstest]
439    fn ct_eq_matches_partial_eq() {
440        let a = Fp::from_u64_reduce(123);
441        let b = Fp::from_u64_reduce(123);
442        let c = Fp::from_u64_reduce(124);
443        assert_eq!(a.ct_eq(b), u64::MAX);
444        assert_eq!(a.ct_eq(c), 0);
445    }
446
447    #[rstest]
448    fn ct_select_picks_branch_by_mask() {
449        let a = Fp::from_u64_reduce(123);
450        let b = Fp::from_u64_reduce(456);
451        assert_eq!(Fp::ct_select(0, a, b), a);
452        assert_eq!(Fp::ct_select(u64::MAX, a, b), b);
453    }
454
455    proptest! {
456        /// Field addition is commutative: `a + b == b + a` for any pair.
457        #[rstest]
458        fn prop_add_commutative(a in arb_fp(), b in arb_fp()) {
459            prop_assert_eq!(a + b, b + a);
460        }
461
462        /// Field addition is associative: `(a + b) + c == a + (b + c)`.
463        #[rstest]
464        fn prop_add_associative(a in arb_fp(), b in arb_fp(), c in arb_fp()) {
465            prop_assert_eq!((a + b) + c, a + (b + c));
466        }
467
468        /// Multiplication distributes over addition.
469        #[rstest]
470        fn prop_distributive(a in arb_fp(), b in arb_fp(), c in arb_fp()) {
471            prop_assert_eq!(a * (b + c), a * b + a * c);
472        }
473
474        /// Multiplication is commutative.
475        #[rstest]
476        fn prop_mul_commutative(a in arb_fp(), b in arb_fp()) {
477            prop_assert_eq!(a * b, b * a);
478        }
479
480        /// Multiplication is associative.
481        #[rstest]
482        fn prop_mul_associative(a in arb_fp(), b in arb_fp(), c in arb_fp()) {
483            prop_assert_eq!((a * b) * c, a * (b * c));
484        }
485
486        /// `a + (-a) == 0` for any element.
487        #[rstest]
488        fn prop_neg_round_trip(a in arb_fp()) {
489            prop_assert_eq!(a + (-a), Fp::ZERO);
490        }
491
492        /// `a - b == a + (-b)` for any pair.
493        #[rstest]
494        fn prop_sub_via_add_neg(a in arb_fp(), b in arb_fp()) {
495            prop_assert_eq!(a - b, a + (-b));
496        }
497
498        /// `(a + b) - b == a` for any pair.
499        #[rstest]
500        fn prop_sub_round_trip(a in arb_fp(), b in arb_fp()) {
501            prop_assert_eq!((a + b) - b, a);
502        }
503
504        /// Squaring matches self-multiplication.
505        #[rstest]
506        fn prop_square_matches_self_mul(a in arb_fp()) {
507            prop_assert_eq!(a.square(), a * a);
508        }
509
510        /// `a * a.invert() == 1` for any non-zero element.
511        #[rstest]
512        fn prop_invert_round_trip(a in arb_fp_nonzero()) {
513            prop_assert_eq!(a * a.invert(), Fp::ONE);
514        }
515
516        /// Fermat's little theorem: `a^(p-1) == 1` for any non-zero element.
517        /// Pins the exponent ladder, the Montgomery reduction, and the
518        /// addition chain through `pow` simultaneously.
519        #[rstest]
520        fn prop_fermat_little(a in arb_fp_nonzero()) {
521            prop_assert_eq!(a.pow(MODULUS - 1), Fp::ONE);
522        }
523
524        /// `(a^2).sqrt()^2 == a^2`: sqrt of a known square round-trips.
525        /// Picks one of the two roots; we only assert the squared identity.
526        #[rstest]
527        fn prop_sqrt_round_trip(a in arb_fp()) {
528            let sq = a.square();
529            let s = sq.sqrt().expect("squares are quadratic residues");
530            prop_assert_eq!(s.square(), sq);
531        }
532
533        /// Canonical bytes round-trip: any `Fp` encodes to canonical bytes
534        /// that decode back to the same element.
535        #[rstest]
536        fn prop_le_bytes_round_trip(a in arb_fp()) {
537            let bytes = a.to_le_bytes();
538            prop_assert_eq!(Fp::try_from_le_bytes(bytes).unwrap(), a);
539        }
540
541        /// `from_u64_canonical` accepts every value in `0..MODULUS`.
542        #[rstest]
543        fn prop_from_u64_canonical_accepts_in_range(v in 0u64..MODULUS) {
544            let f = Fp::from_u64_canonical(v).expect("in-range value");
545            prop_assert_eq!(f.to_u64(), v);
546        }
547
548        /// `from_u64_canonical` rejects every value `>= MODULUS`.
549        #[rstest]
550        fn prop_from_u64_canonical_rejects_out_of_range(v in MODULUS..=u64::MAX) {
551            prop_assert!(Fp::from_u64_canonical(v).is_none());
552        }
553
554        /// `msquare(n)` matches `n` iterated `square` calls.
555        #[rstest]
556        fn prop_msquare_matches_iterated_square(a in arb_fp(), n in 0u32..16) {
557            let mut iter = a;
558            for _ in 0..n {
559                iter = iter.square();
560            }
561            prop_assert_eq!(a.msquare(n), iter);
562        }
563
564        /// `ct_eq` agrees with `==` over arbitrary pairs.
565        #[rstest]
566        fn prop_ct_eq_matches_partial_eq(a in arb_fp(), b in arb_fp()) {
567            let ct = a.ct_eq(b);
568            if a == b {
569                prop_assert_eq!(ct, u64::MAX);
570            } else {
571                prop_assert_eq!(ct, 0);
572            }
573        }
574
575        /// `ct_select` picks `a` for mask 0 and `b` for mask u64::MAX.
576        #[rstest]
577        fn prop_ct_select_picks_branch(a in arb_fp(), b in arb_fp()) {
578            prop_assert_eq!(Fp::ct_select(0, a, b), a);
579            prop_assert_eq!(Fp::ct_select(u64::MAX, a, b), b);
580        }
581    }
582
583    #[rstest]
584    fn matches_go_reference_vectors() {
585        let suite: Vectors = serde_json::from_str(VECTORS_JSON).expect("parse vectors");
586        assert!(!suite.vectors.is_empty(), "vector file is empty");
587
588        for (i, v) in suite.vectors.iter().enumerate() {
589            let a = Fp::try_from_le_bytes(decode_le8(&v.a))
590                .unwrap_or_else(|| panic!("vector {i}: decode a"));
591            let b = Fp::try_from_le_bytes(decode_le8(&v.b))
592                .unwrap_or_else(|| panic!("vector {i}: decode b"));
593            let e = parse_u64(&v.e);
594
595            assert_eq!((a + b).to_le_bytes(), decode_le8(&v.add), "vector {i}: add");
596            assert_eq!((a - b).to_le_bytes(), decode_le8(&v.sub), "vector {i}: sub");
597            assert_eq!((a * b).to_le_bytes(), decode_le8(&v.mul), "vector {i}: mul");
598            assert_eq!((-a).to_le_bytes(), decode_le8(&v.neg_a), "vector {i}: neg");
599            assert_eq!(
600                a.invert().to_le_bytes(),
601                decode_le8(&v.inv_a),
602                "vector {i}: inv"
603            );
604            assert_eq!(
605                a.pow(e).to_le_bytes(),
606                decode_le8(&v.pow_a_e),
607                "vector {i}: pow"
608            );
609            assert_eq!(a == b, v.a_eq_b, "vector {i}: eq");
610        }
611    }
612}
nautilus_lighter/signing/field/goldilocks.rs

nautilus_lighter/signing/field/
goldilocks.rs
