Goldilocks field (#227)

* Goldilocks field

Based on Hamish's old branch, but I updated it with a few missing things like generators.

Pulled the inversion code into a shared helper method to avoid redundancy.

Just the base field for now. We can add a quartic extension field later.

* typo

* PR feedback

* More overflowing -> wrapping

* fmt

* cleanup

benches/field_arithmetic.rs

@@ -4,6 +4,7 @@ use criterion::{criterion_group, criterion_main, BatchSize, Criterion};
 use plonky2::field::crandall_field::CrandallField;
 use plonky2::field::extension_field::quartic::QuarticCrandallField;
 use plonky2::field::field_types::Field;
+use plonky2::field::goldilocks_field::GoldilocksField;
 use tynm::type_name;
 
 pub(crate) fn bench_field<F: Field>(c: &mut Criterion) {
@@ -90,6 +91,7 @@ pub(crate) fn bench_field<F: Field>(c: &mut Criterion) {
 
 fn criterion_benchmark(c: &mut Criterion) {
     bench_field::<CrandallField>(c);
+    bench_field::<GoldilocksField>(c);
     bench_field::<QuarticCrandallField>(c);
 }
 

src/field/crandall_field.rs

@@ -5,7 +5,6 @@ use std::iter::{Product, Sum};
 use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
 use num::bigint::BigUint;
-use num::Integer;
 use rand::Rng;
 use serde::{Deserialize, Serialize};
 
@@ -13,6 +12,7 @@ use crate::field::extension_field::quadratic::QuadraticCrandallField;
 use crate::field::extension_field::quartic::QuarticCrandallField;
 use crate::field::extension_field::{Extendable, Frobenius};
 use crate::field::field_types::{Field, PrimeField, RichField};
+use crate::field::inversion::try_inverse_u64;
 
 /// EPSILON = 9 * 2**28 - 1
 const EPSILON: u64 = 2415919103;
@@ -160,82 +160,11 @@ impl Field for CrandallField {
     const POWER_OF_TWO_GENERATOR: Self = Self(10281950781551402419);
 
     fn order() -> BigUint {
-        BigUint::from(Self::ORDER)
+        Self::ORDER.into()
     }
 
-    #[inline]
-    fn square(&self) -> Self {
-        *self * *self
-    }
-
-    #[inline]
-    fn cube(&self) -> Self {
-        *self * *self * *self
-    }
-
-    #[allow(clippy::many_single_char_names)] // The names are from the paper.
     fn try_inverse(&self) -> Option<Self> {
-        if self.is_zero() {
-            return None;
-        }
-
-        // Based on Algorithm 16 of "Efficient Software-Implementation of Finite Fields with
-        // Applications to Cryptography".
-
-        let p = Self::ORDER;
-        let mut u = self.to_canonical_u64();
-        let mut v = p;
-        let mut b = 1u64;
-        let mut c = 0u64;
-
-        while u != 1 && v != 1 {
-            let u_tz = u.trailing_zeros();
-            u >>= u_tz;
-            for _ in 0..u_tz {
-                if b.is_even() {
-                    b /= 2;
-                } else {
-                    // b = (b + p)/2, avoiding overflow
-                    b = (b / 2) + (p / 2) + 1;
-                }
-            }
-
-            let v_tz = v.trailing_zeros();
-            v >>= v_tz;
-            for _ in 0..v_tz {
-                if c.is_even() {
-                    c /= 2;
-                } else {
-                    // c = (c + p)/2, avoiding overflow
-                    c = (c / 2) + (p / 2) + 1;
-                }
-            }
-
-            if u >= v {
-                u -= v;
-                // b -= c
-                let (mut diff, under) = b.overflowing_sub(c);
-                if under {
-                    diff = diff.overflowing_add(p).0;
-                }
-                b = diff;
-            } else {
-                v -= u;
-                // c -= b
-                let (mut diff, under) = c.overflowing_sub(b);
-                if under {
-                    diff = diff.overflowing_add(p).0;
-                }
-                c = diff;
-            }
-        }
-
-        let inverse = Self(if u == 1 { b } else { c });
-
-        // Should change to debug_assert_eq; using assert_eq as an extra precaution for now until
-        // we're more confident the impl is correct.
-        assert_eq!(*self * inverse, Self::ONE);
-        Some(inverse)
+        try_inverse_u64(self.0, Self::ORDER).map(|inv| Self(inv))
     }
 
     #[inline]
@@ -380,7 +309,7 @@ impl Add for CrandallField {
     #[allow(clippy::suspicious_arithmetic_impl)]
     fn add(self, rhs: Self) -> Self {
         let (sum, over) = self.0.overflowing_add(rhs.to_canonical_u64());
-        Self(sum.overflowing_sub((over as u64) * Self::ORDER).0)
+        Self(sum.wrapping_sub((over as u64) * Self::ORDER))
     }
 }
 
@@ -403,7 +332,7 @@ impl Sub for CrandallField {
     #[allow(clippy::suspicious_arithmetic_impl)]
     fn sub(self, rhs: Self) -> Self {
         let (diff, under) = self.0.overflowing_sub(rhs.to_canonical_u64());
-        Self(diff.overflowing_add((under as u64) * Self::ORDER).0)
+        Self(diff.wrapping_add((under as u64) * Self::ORDER))
     }
 }
 
@@ -470,7 +399,7 @@ impl RichField for CrandallField {}
 #[inline]
 unsafe fn add_no_canonicalize(lhs: CrandallField, rhs: CrandallField) -> CrandallField {
     let (sum, over) = lhs.0.overflowing_add(rhs.0);
-    CrandallField(sum.overflowing_sub((over as u64) * CrandallField::ORDER).0)
+    CrandallField(sum.wrapping_sub((over as u64) * CrandallField::ORDER))
 }
 
 /// Reduces to a 64-bit value. The result might not be in canonical form; it could be in between the

src/field/goldilocks_field.rs

@@ -0,0 +1,240 @@
+use std::fmt;
+use std::fmt::{Debug, Display, Formatter};
+use std::hash::{Hash, Hasher};
+use std::iter::{Product, Sum};
+use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
+
+use num::BigUint;
+use rand::Rng;
+use serde::{Deserialize, Serialize};
+
+use crate::field::field_types::{Field, PrimeField};
+use crate::field::inversion::try_inverse_u64;
+
+const EPSILON: u64 = (1 << 32) - 1;
+
+/// A field selected to have fast reduction.
+///
+/// Its order is 2^64 - 2^32 + 1.
+/// ```ignore
+/// P = 2**64 - EPSILON
+///   = 2**64 - 2**32 + 1
+///   = 2**32 * (2**32 - 1) + 1
+/// ```
+#[derive(Copy, Clone, Serialize, Deserialize)]
+pub struct GoldilocksField(pub u64);
+
+impl Default for GoldilocksField {
+    fn default() -> Self {
+        Self::ZERO
+    }
+}
+
+impl PartialEq for GoldilocksField {
+    fn eq(&self, other: &Self) -> bool {
+        self.to_canonical_u64() == other.to_canonical_u64()
+    }
+}
+
+impl Eq for GoldilocksField {}
+
+impl Hash for GoldilocksField {
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        state.write_u64(self.to_canonical_u64())
+    }
+}
+
+impl Display for GoldilocksField {
+    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
+        Display::fmt(&self.0, f)
+    }
+}
+
+impl Debug for GoldilocksField {
+    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
+        Debug::fmt(&self.0, f)
+    }
+}
+
+impl Field for GoldilocksField {
+    type PrimeField = Self;
+
+    const ZERO: Self = Self(0);
+    const ONE: Self = Self(1);
+    const TWO: Self = Self(2);
+    const NEG_ONE: Self = Self(Self::ORDER - 1);
+    const CHARACTERISTIC: u64 = Self::ORDER;
+
+    const TWO_ADICITY: usize = 32;
+
+    // Sage: `g = GF(p).multiplicative_generator()`
+    const MULTIPLICATIVE_GROUP_GENERATOR: Self = Self(7);
+
+    // Sage:
+    // ```
+    // g_2 = g^((p - 1) / 2^32)
+    // g_2.multiplicative_order()
+    // ```
+    const POWER_OF_TWO_GENERATOR: Self = Self(1753635133440165772);
+
+    fn order() -> BigUint {
+        Self::ORDER.into()
+    }
+
+    fn try_inverse(&self) -> Option<Self> {
+        try_inverse_u64(self.0, Self::ORDER).map(|inv| Self(inv))
+    }
+
+    #[inline]
+    fn from_canonical_u64(n: u64) -> Self {
+        Self(n)
+    }
+
+    fn from_noncanonical_u128(n: u128) -> Self {
+        reduce128(n)
+    }
+
+    fn rand_from_rng<R: Rng>(rng: &mut R) -> Self {
+        Self::from_canonical_u64(rng.gen_range(0..Self::ORDER))
+    }
+}
+
+impl PrimeField for GoldilocksField {
+    const ORDER: u64 = 0xFFFFFFFF00000001;
+
+    #[inline]
+    fn to_canonical_u64(&self) -> u64 {
+        let mut c = self.0;
+        // We only need one condition subtraction, since 2 * ORDER would not fit in a u64.
+        if c >= Self::ORDER {
+            c -= Self::ORDER;
+        }
+        c
+    }
+
+    fn to_noncanonical_u64(&self) -> u64 {
+        self.0
+    }
+}
+
+impl Neg for GoldilocksField {
+    type Output = Self;
+
+    #[inline]
+    fn neg(self) -> Self {
+        if self.is_zero() {
+            Self::ZERO
+        } else {
+            Self(Self::ORDER - self.to_canonical_u64())
+        }
+    }
+}
+
+impl Add for GoldilocksField {
+    type Output = Self;
+
+    #[inline]
+    #[allow(clippy::suspicious_arithmetic_impl)]
+    fn add(self, rhs: Self) -> Self {
+        let (sum, over) = self.0.overflowing_add(rhs.to_canonical_u64());
+        Self(sum.wrapping_sub((over as u64) * Self::ORDER))
+    }
+}
+
+impl AddAssign for GoldilocksField {
+    fn add_assign(&mut self, rhs: Self) {
+        *self = *self + rhs;
+    }
+}
+
+impl Sum for GoldilocksField {
+    fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
+        iter.fold(Self::ZERO, |acc, x| acc + x)
+    }
+}
+
+impl Sub for GoldilocksField {
+    type Output = Self;
+
+    #[inline]
+    #[allow(clippy::suspicious_arithmetic_impl)]
+    fn sub(self, rhs: Self) -> Self {
+        let (diff, under) = self.0.overflowing_sub(rhs.to_canonical_u64());
+        Self(diff.wrapping_add((under as u64) * Self::ORDER))
+    }
+}
+
+impl SubAssign for GoldilocksField {
+    #[inline]
+    fn sub_assign(&mut self, rhs: Self) {
+        *self = *self - rhs;
+    }
+}
+
+impl Mul for GoldilocksField {
+    type Output = Self;
+
+    #[inline]
+    fn mul(self, rhs: Self) -> Self {
+        reduce128((self.0 as u128) * (rhs.0 as u128))
+    }
+}
+
+impl MulAssign for GoldilocksField {
+    #[inline]
+    fn mul_assign(&mut self, rhs: Self) {
+        *self = *self * rhs;
+    }
+}
+
+impl Product for GoldilocksField {
+    fn product<I: Iterator<Item = Self>>(iter: I) -> Self {
+        iter.fold(Self::ONE, |acc, x| acc * x)
+    }
+}
+
+impl Div for GoldilocksField {
+    type Output = Self;
+
+    #[allow(clippy::suspicious_arithmetic_impl)]
+    fn div(self, rhs: Self) -> Self::Output {
+        self * rhs.inverse()
+    }
+}
+
+impl DivAssign for GoldilocksField {
+    fn div_assign(&mut self, rhs: Self) {
+        *self = *self / rhs;
+    }
+}
+
+/// Reduces to a 64-bit value. The result might not be in canonical form; it could be in between the
+/// field order and `2^64`.
+#[inline]
+fn reduce128(x: u128) -> GoldilocksField {
+    let (x_lo, x_hi) = split(x); // This is a no-op
+    let x_hi_hi = x_hi >> 32;
+    let x_hi_lo = x_hi & EPSILON;
+
+    let (mut t0, borrow) = x_lo.overflowing_sub(x_hi_hi);
+    t0 = t0.wrapping_sub(EPSILON * (borrow as u64));
+
+    let t1 = x_hi_lo * EPSILON;
+
+    let (mut t2, carry) = t1.overflowing_add(t0);
+    t2 = t2.wrapping_add(EPSILON * (carry as u64));
+    GoldilocksField(t2)
+}
+
+#[inline]
+fn split(x: u128) -> (u64, u64) {
+    (x as u64, (x >> 64) as u64)
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::{test_field_arithmetic, test_prime_field_arithmetic};
+
+    test_prime_field_arithmetic!(crate::field::goldilocks_field::GoldilocksField);
+    test_field_arithmetic!(crate::field::goldilocks_field::GoldilocksField);
+}

src/field/inversion.rs

@@ -0,0 +1,61 @@
+use num::{Integer, Zero};
+
+/// Try to invert an element in a prime field with the given modulus.
+#[allow(clippy::many_single_char_names)] // The names are from the paper.
+pub(crate) fn try_inverse_u64(x: u64, p: u64) -> Option<u64> {
+    if x.is_zero() {
+        return None;
+    }
+
+    // Based on Algorithm 16 of "Efficient Software-Implementation of Finite Fields with
+    // Applications to Cryptography".
+
+    let mut u = x;
+    let mut v = p;
+    let mut b = 1u64;
+    let mut c = 0u64;
+
+    while u != 1 && v != 1 {
+        let u_tz = u.trailing_zeros();
+        u >>= u_tz;
+        for _ in 0..u_tz {
+            if b.is_even() {
+                b /= 2;
+            } else {
+                // b = (b + p)/2, avoiding overflow
+                b = (b / 2) + (p / 2) + 1;
+            }
+        }
+
+        let v_tz = v.trailing_zeros();
+        v >>= v_tz;
+        for _ in 0..v_tz {
+            if c.is_even() {
+                c /= 2;
+            } else {
+                // c = (c + p)/2, avoiding overflow
+                c = (c / 2) + (p / 2) + 1;
+            }
+        }
+
+        if u >= v {
+            u -= v;
+            // b -= c
+            let (mut diff, under) = b.overflowing_sub(c);
+            if under {
+                diff = diff.wrapping_add(p);
+            }
+            b = diff;
+        } else {
+            v -= u;
+            // c -= b
+            let (mut diff, under) = c.overflowing_sub(b);
+            if under {
+                diff = diff.wrapping_add(p);
+            }
+            c = diff;
+        }
+    }
+
+    Some(if u == 1 { b } else { c })
+}

src/field/mod.rs

@@ -3,7 +3,9 @@ pub mod crandall_field;
 pub mod extension_field;
 pub mod fft;
 pub mod field_types;
+pub mod goldilocks_field;
 pub(crate) mod interpolation;
+mod inversion;
 pub(crate) mod packable;
 pub(crate) mod packed_field;
 

src/field/packed_crandall_avx2.rs

@@ -198,7 +198,7 @@ impl Sum for PackedCrandallAVX2 {
 }
 
 const EPSILON: u64 = (1 << 31) + (1 << 28) - 1;
-const FIELD_ORDER: u64 = 0u64.overflowing_sub(EPSILON).0;
+const FIELD_ORDER: u64 = 0u64.wrapping_sub(EPSILON);
 const SIGN_BIT: u64 = 1 << 63;
 
 #[inline]