Quintic extension fields (#489)

* Initial implementation of quintic extensions.

* Update to/from_biguint() methods.

* cargo fmt

* Fix call to test suite.

* Small optimisation in try_inverse().

* Replace multiplicative group generator and document requirement.

field/src/extension_field/mod.rs

@@ -5,6 +5,7 @@ use crate::field_types::Field;
 pub mod algebra;
 pub mod quadratic;
 pub mod quartic;
+pub mod quintic;
 
 /// Optimal extension field trait.
 /// A degree `d` field extension is optimal if there exists a base field element `W`,
@@ -67,6 +68,8 @@ pub trait Extendable<const D: usize>: Field + Sized {
 
     const DTH_ROOT: Self;
 
+    /// Chosen so that when raised to the power `(p^D - 1) >> F::Extension::TWO_ADICITY)`
+    /// we obtain F::EXT_POWER_OF_TWO_GENERATOR.
     const EXT_MULTIPLICATIVE_GROUP_GENERATOR: [Self; D];
 
     /// Chosen so that when raised to the power `1<<(Self::TWO_ADICITY-Self::BaseField::TWO_ADICITY)`,

field/src/extension_field/quintic.rs

@@ -0,0 +1,278 @@
+use std::fmt::{Debug, Display, Formatter};
+use std::iter::{Product, Sum};
+use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
+
+use num::bigint::BigUint;
+use num::traits::Pow;
+use rand::Rng;
+use serde::{Deserialize, Serialize};
+
+use crate::extension_field::{Extendable, FieldExtension, Frobenius, OEF};
+use crate::field_types::Field;
+use crate::ops::Square;
+
+#[derive(Copy, Clone, Eq, PartialEq, Hash, Serialize, Deserialize)]
+#[serde(bound = "")]
+pub struct QuinticExtension<F: Extendable<5>>(pub [F; 5]);
+
+impl<F: Extendable<5>> Default for QuinticExtension<F> {
+    fn default() -> Self {
+        Self::ZERO
+    }
+}
+
+impl<F: Extendable<5>> OEF<5> for QuinticExtension<F> {
+    const W: F = F::W;
+    const DTH_ROOT: F = F::DTH_ROOT;
+}
+
+impl<F: Extendable<5>> Frobenius<5> for QuinticExtension<F> {}
+
+impl<F: Extendable<5>> FieldExtension<5> for QuinticExtension<F> {
+    type BaseField = F;
+
+    fn to_basefield_array(&self) -> [F; 5] {
+        self.0
+    }
+
+    fn from_basefield_array(arr: [F; 5]) -> Self {
+        Self(arr)
+    }
+
+    fn from_basefield(x: F) -> Self {
+        x.into()
+    }
+}
+
+impl<F: Extendable<5>> From<F> for QuinticExtension<F> {
+    fn from(x: F) -> Self {
+        Self([x, F::ZERO, F::ZERO, F::ZERO, F::ZERO])
+    }
+}
+
+impl<F: Extendable<5>> Field for QuinticExtension<F> {
+    const ZERO: Self = Self([F::ZERO; 5]);
+    const ONE: Self = Self([F::ONE, F::ZERO, F::ZERO, F::ZERO, F::ZERO]);
+    const TWO: Self = Self([F::TWO, F::ZERO, F::ZERO, F::ZERO, F::ZERO]);
+    const NEG_ONE: Self = Self([F::NEG_ONE, F::ZERO, F::ZERO, F::ZERO, F::ZERO]);
+
+    // `p^5 - 1 = (p - 1)(p^4 + p^3 + p^2 + p + 1)`. The `p - 1` term
+    // has a two-adicity of `F::TWO_ADICITY` and the term `p^4 + p^3 +
+    // p^2 + p + 1` is odd since it is the sum of an odd number of odd
+    // terms. Hence the two-adicity of `p^5 - 1` is the same as for
+    // `p - 1`.
+    const TWO_ADICITY: usize = F::TWO_ADICITY;
+    const CHARACTERISTIC_TWO_ADICITY: usize = F::CHARACTERISTIC_TWO_ADICITY;
+
+    const MULTIPLICATIVE_GROUP_GENERATOR: Self = Self(F::EXT_MULTIPLICATIVE_GROUP_GENERATOR);
+    const POWER_OF_TWO_GENERATOR: Self = Self(F::EXT_POWER_OF_TWO_GENERATOR);
+
+    const BITS: usize = F::BITS * 5;
+
+    fn order() -> BigUint {
+        F::order().pow(5u32)
+    }
+    fn characteristic() -> BigUint {
+        F::characteristic()
+    }
+
+    // Algorithm 11.3.4 in Handbook of Elliptic and Hyperelliptic Curve Cryptography.
+    fn try_inverse(&self) -> Option<Self> {
+        if self.is_zero() {
+            return None;
+        }
+
+        // Writing 'a' for self:
+        let d = self.frobenius(); // d = a^p
+        let e = d * d.frobenius(); // e = a^(p + p^2)
+        let f = e * e.repeated_frobenius(2); // f = a^(p + p^2 + p^3 + p^4)
+
+        // f contains a^(r-1) and a^r is in the base field.
+        debug_assert!(FieldExtension::<5>::is_in_basefield(&(*self * f)));
+
+        // g = a^r is in the base field, so only compute that
+        // coefficient rather than the full product. The equation is
+        // extracted from Mul::mul(...) below.
+        let Self([a0, a1, a2, a3, a4]) = *self;
+        let Self([b0, b1, b2, b3, b4]) = f;
+        let g = a0 * b0 + <Self as OEF<5>>::W * (a1 * b4 + a2 * b3 + a3 * b2 + a4 * b1);
+
+        Some(FieldExtension::<5>::scalar_mul(&f, g.inverse()))
+    }
+
+    fn from_biguint(n: BigUint) -> Self {
+        Self([F::from_biguint(n), F::ZERO, F::ZERO, F::ZERO, F::ZERO])
+    }
+
+    fn from_canonical_u64(n: u64) -> Self {
+        F::from_canonical_u64(n).into()
+    }
+
+    fn from_noncanonical_u128(n: u128) -> Self {
+        F::from_noncanonical_u128(n).into()
+    }
+
+    fn rand_from_rng<R: Rng>(rng: &mut R) -> Self {
+        Self::from_basefield_array([
+            F::rand_from_rng(rng),
+            F::rand_from_rng(rng),
+            F::rand_from_rng(rng),
+            F::rand_from_rng(rng),
+            F::rand_from_rng(rng),
+        ])
+    }
+}
+
+impl<F: Extendable<5>> Display for QuinticExtension<F> {
+    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
+        write!(
+            f,
+            "{} + {}*a + {}*a^2 + {}*a^3 + {}*a^4",
+            self.0[0], self.0[1], self.0[2], self.0[3], self.0[4]
+        )
+    }
+}
+
+impl<F: Extendable<5>> Debug for QuinticExtension<F> {
+    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
+        Display::fmt(self, f)
+    }
+}
+
+impl<F: Extendable<5>> Neg for QuinticExtension<F> {
+    type Output = Self;
+
+    #[inline]
+    fn neg(self) -> Self {
+        Self([-self.0[0], -self.0[1], -self.0[2], -self.0[3], -self.0[4]])
+    }
+}
+
+impl<F: Extendable<5>> Add for QuinticExtension<F> {
+    type Output = Self;
+
+    #[inline]
+    fn add(self, rhs: Self) -> Self {
+        Self([
+            self.0[0] + rhs.0[0],
+            self.0[1] + rhs.0[1],
+            self.0[2] + rhs.0[2],
+            self.0[3] + rhs.0[3],
+            self.0[4] + rhs.0[4],
+        ])
+    }
+}
+
+impl<F: Extendable<5>> AddAssign for QuinticExtension<F> {
+    fn add_assign(&mut self, rhs: Self) {
+        *self = *self + rhs;
+    }
+}
+
+impl<F: Extendable<5>> Sum for QuinticExtension<F> {
+    fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
+        iter.fold(Self::ZERO, |acc, x| acc + x)
+    }
+}
+
+impl<F: Extendable<5>> Sub for QuinticExtension<F> {
+    type Output = Self;
+
+    #[inline]
+    fn sub(self, rhs: Self) -> Self {
+        Self([
+            self.0[0] - rhs.0[0],
+            self.0[1] - rhs.0[1],
+            self.0[2] - rhs.0[2],
+            self.0[3] - rhs.0[3],
+            self.0[4] - rhs.0[4],
+        ])
+    }
+}
+
+impl<F: Extendable<5>> SubAssign for QuinticExtension<F> {
+    #[inline]
+    fn sub_assign(&mut self, rhs: Self) {
+        *self = *self - rhs;
+    }
+}
+
+impl<F: Extendable<5>> Mul for QuinticExtension<F> {
+    type Output = Self;
+
+    #[inline]
+    fn mul(self, rhs: Self) -> Self {
+        let Self([a0, a1, a2, a3, a4]) = self;
+        let Self([b0, b1, b2, b3, b4]) = rhs;
+        let w = <Self as OEF<5>>::W;
+
+        let c0 = a0 * b0 + w * (a1 * b4 + a2 * b3 + a3 * b2 + a4 * b1);
+        let c1 = a0 * b1 + a1 * b0 + w * (a2 * b4 + a3 * b3 + a4 * b2);
+        let c2 = a0 * b2 + a1 * b1 + a2 * b0 + w * (a3 * b4 + a4 * b3);
+        let c3 = a0 * b3 + a1 * b2 + a2 * b1 + a3 * b0 + w * a4 * b4;
+        let c4 = a0 * b4 + a1 * b3 + a2 * b2 + a3 * b1 + a4 * b0;
+
+        Self([c0, c1, c2, c3, c4])
+    }
+}
+
+impl<F: Extendable<5>> MulAssign for QuinticExtension<F> {
+    #[inline]
+    fn mul_assign(&mut self, rhs: Self) {
+        *self = *self * rhs;
+    }
+}
+
+impl<F: Extendable<5>> Square for QuinticExtension<F> {
+    #[inline(always)]
+    fn square(&self) -> Self {
+        let Self([a0, a1, a2, a3, a4]) = *self;
+        let w = <Self as OEF<5>>::W;
+        let double_w = <Self as OEF<5>>::W.double();
+
+        let c0 = a0.square() + double_w * (a1 * a4 + a2 * a3);
+        let double_a0 = a0.double();
+        let c1 = double_a0 * a1 + double_w * a2 * a4 + w * a3 * a3;
+        let c2 = double_a0 * a2 + a1 * a1 + double_w * a4 * a3;
+        let double_a1 = a1.double();
+        let c3 = double_a0 * a3 + double_a1 * a2 + w * a4 * a4;
+        let c4 = double_a0 * a4 + double_a1 * a3 + a2 * a2;
+
+        Self([c0, c1, c2, c3, c4])
+    }
+}
+
+impl<F: Extendable<5>> Product for QuinticExtension<F> {
+    fn product<I: Iterator<Item = Self>>(iter: I) -> Self {
+        iter.fold(Self::ONE, |acc, x| acc * x)
+    }
+}
+
+impl<F: Extendable<5>> Div for QuinticExtension<F> {
+    type Output = Self;
+
+    #[allow(clippy::suspicious_arithmetic_impl)]
+    fn div(self, rhs: Self) -> Self::Output {
+        self * rhs.inverse()
+    }
+}
+
+impl<F: Extendable<5>> DivAssign for QuinticExtension<F> {
+    fn div_assign(&mut self, rhs: Self) {
+        *self = *self / rhs;
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    mod goldilocks {
+        use crate::{test_field_arithmetic, test_field_extension};
+
+        test_field_extension!(crate::goldilocks_field::GoldilocksField, 5);
+        test_field_arithmetic!(
+            crate::extension_field::quintic::QuinticExtension<
+                crate::goldilocks_field::GoldilocksField,
+            >
+        );
+    }
+}

field/src/goldilocks_field.rs

@@ -11,6 +11,7 @@ use serde::{Deserialize, Serialize};
 
 use crate::extension_field::quadratic::QuadraticExtension;
 use crate::extension_field::quartic::QuarticExtension;
+use crate::extension_field::quintic::QuinticExtension;
 use crate::extension_field::{Extendable, Frobenius};
 use crate::field_types::{Field, Field64, PrimeField, PrimeField64};
 use crate::inversion::try_inverse_u64;
@@ -317,6 +318,31 @@ impl Extendable<4> for GoldilocksField {
         [Self(0), Self(0), Self(0), Self(12587610116473453104)];
 }
 
+impl Extendable<5> for GoldilocksField {
+    type Extension = QuinticExtension<Self>;
+
+    const W: Self = Self(3);
+
+    // DTH_ROOT = W^((ORDER - 1)/5)
+    const DTH_ROOT: Self = Self(1041288259238279555);
+
+    const EXT_MULTIPLICATIVE_GROUP_GENERATOR: [Self; 5] = [
+        Self(2899034827742553394),
+        Self(13012057356839176729),
+        Self(14593811582388663055),
+        Self(7722900811313895436),
+        Self(4557222484695340057),
+    ];
+
+    const EXT_POWER_OF_TWO_GENERATOR: [Self; 5] = [
+        Self::POWER_OF_TWO_GENERATOR,
+        Self(0),
+        Self(0),
+        Self(0),
+        Self(0),
+    ];
+}
+
 /// Fast addition modulo ORDER for x86-64.
 /// This function is marked unsafe for the following reasons:
 ///   - It is only correct if x + y < 2**64 + ORDER = 0x1ffffffff00000001.

plonky2/benches/field_arithmetic.rs

@@ -1,5 +1,6 @@
 use criterion::{criterion_group, criterion_main, BatchSize, Criterion};
 use plonky2::field::extension_field::quartic::QuarticExtension;
+use plonky2::field::extension_field::quintic::QuinticExtension;
 use plonky2::field::field_types::Field;
 use plonky2::field::goldilocks_field::GoldilocksField;
 use tynm::type_name;
@@ -175,6 +176,7 @@ pub(crate) fn bench_field<F: Field>(c: &mut Criterion) {
 fn criterion_benchmark(c: &mut Criterion) {
     bench_field::<GoldilocksField>(c);
     bench_field::<QuarticExtension<GoldilocksField>>(c);
+    bench_field::<QuinticExtension<GoldilocksField>>(c);
 }
 
 criterion_group!(benches, criterion_benchmark);