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

nautilus_lighter/signing/field/
goldilocks.rs
