diff --git a/examples/dkg.rs b/examples/dkg.rs index ed80e0b..737ce07 100644 --- a/examples/dkg.rs +++ b/examples/dkg.rs @@ -165,7 +165,7 @@ fn main() -> Result<(), confy::ConfyError> { tracing_subscriber::fmt().init(); event!(Level::INFO, "Begin dkg::main"); let config_path = PathBuf::from("examples/dkg.toml"); - let cfg: MyConfig = confy::load_path(&config_path)?; + let cfg: MyConfig = confy::load_path(config_path)?; event!(Level::INFO, "Loaded config: {:?}", cfg); const NUMBER_OF_ROUNDS: u64 = 3; diff --git a/src/fields/extensions.rs b/src/fields/extensions.rs index 22d2649..1c53604 100644 --- a/src/fields/extensions.rs +++ b/src/fields/extensions.rs @@ -72,7 +72,7 @@ impl<'a, 'b, const D: usize, const N: usize, F: FieldExtensionTrait> Add<&'b FieldExtension> for &'a FieldExtension { type Output = FieldExtension; - + #[inline] fn add(self, other: &'b FieldExtension) -> Self::Output { let mut i = 0; let mut retval = [F::zero(); N]; @@ -87,6 +87,7 @@ impl> Add { type Output = Self; + #[inline] fn add(self, other: FieldExtension) -> Self::Output { &self + &other } @@ -94,6 +95,7 @@ impl> Add> AddAssign for FieldExtension { + #[inline] fn add_assign(&mut self, other: Self) { *self = *self + other; } @@ -102,7 +104,7 @@ impl<'a, 'b, const D: usize, const N: usize, F: FieldExtensionTrait> Sub<&'b FieldExtension> for &'a FieldExtension { type Output = FieldExtension; - + #[inline] fn sub(self, other: &'b FieldExtension) -> Self::Output { let mut i = 0; let mut retval = [F::zero(); N]; @@ -117,6 +119,7 @@ impl> Sub { type Output = Self; + #[inline] fn sub(self, other: FieldExtension) -> Self::Output { &self - &other } @@ -124,6 +127,7 @@ impl> Sub> SubAssign for FieldExtension { + #[inline] fn sub_assign(&mut self, other: Self) { *self = *self - other; } @@ -145,6 +149,7 @@ impl> PartialEq } impl> Neg for FieldExtension { type Output = Self; + #[inline] fn neg(self) -> Self { let mut i = 0; let mut retval = [F::zero(); N]; @@ -158,6 +163,7 @@ impl> Neg for Field impl> Zero for FieldExtension { + #[inline] fn zero() -> Self { Self::new(&[F::zero(); N]) } diff --git a/src/fields/fp.rs b/src/fields/fp.rs index bfd0794..b834657 100644 --- a/src/fields/fp.rs +++ b/src/fields/fp.rs @@ -49,6 +49,12 @@ pub(crate) const BN254_FP_MODULUS: Fp = Fp::new(U256::from_words([ 0xB85045B68181585D, 0x30644E72E131A029, ])); +pub(crate) const FP_QUADRATIC_NON_RESIDUE: Fp = Fp::new(U256::from_words([ + 4332616871279656262, + 10917124144477883021, + 13281191951274694749, + 3486998266802970665, +])); /// This defines the key properties of a field extension. Now, mathematically, /// a finite field satisfies many rigorous mathematical properties. The /// (non-exhaustive) list below simply suffices to illustrate those properties @@ -83,9 +89,6 @@ pub trait FieldExtensionTrait: + Inv + From { - /// multiplication in a field extension is dictated heavily such a value below, - /// so we define the quadratic non-residue here - fn quadratic_non_residue() -> Self; /// generate a random value in the field extension based on the random number generator from /// `crypto_bigint` fn rand(rng: &mut R) -> Self; @@ -200,11 +203,6 @@ macro_rules! define_finite_prime_field { // appropriate degree, in our case degree 1 (with // therefore 1 unique representation of an element) impl FieldExtensionTrait<$degree, $nreps> for $wrapper_name { - fn quadratic_non_residue() -> Self { - //this is p - 1 mod p = -1 mod p = 0 - 1 mod p - // = -1 - Self::new((-Self::ONE).1.retrieve()) - } /// Generate a random value in the field /// # Arguments /// * `rng` - R: CryptoRngCore - the random number generator to use @@ -233,11 +231,13 @@ macro_rules! define_finite_prime_field { /// of the field element. All binops with assignment equivalents are given impl Add for $wrapper_name { type Output = Self; + #[inline] fn add(self, other: Self) -> Self { Self::new((self.1 + other.1).retrieve()) } } impl AddAssign for $wrapper_name { + #[inline] fn add_assign(&mut self, other: Self) { *self = *self + other; } @@ -262,11 +262,14 @@ macro_rules! define_finite_prime_field { } impl Sub for $wrapper_name { type Output = Self; + + #[inline] fn sub(self, other: Self) -> Self { Self::new((self.1 - other.1).retrieve()) } } impl SubAssign for $wrapper_name { + #[inline] fn sub_assign(&mut self, other: Self) { *self = *self - other; } @@ -301,11 +304,13 @@ macro_rules! define_finite_prime_field { } impl Mul for $wrapper_name { type Output = Self; + #[inline] fn mul(self, other: Self) -> Self { Self::new((self.1 * other.1).retrieve()) } } impl MulAssign for $wrapper_name { + #[inline] fn mul_assign(&mut self, other: Self) { *self = *self * other; } @@ -324,6 +329,7 @@ macro_rules! define_finite_prime_field { /// impl Inv for $wrapper_name { type Output = Self; + #[inline] fn inv(self) -> Self { Self::new((CtOption::from(self.1.inv()).unwrap_or(Self::from(0u64).1)).retrieve()) } @@ -331,23 +337,28 @@ macro_rules! define_finite_prime_field { #[allow(clippy::suspicious_arithmetic_impl)] impl Div for $wrapper_name { type Output = Self; + #[inline] fn div(self, other: Self) -> Self { self * other.inv() } } impl DivAssign for $wrapper_name { + #[inline] fn div_assign(&mut self, other: Self) { *self = *self / other; } } impl Neg for $wrapper_name { type Output = Self; + + #[inline] fn neg(self) -> Self { Self::new((-self.1).retrieve()) } } impl Pow for $wrapper_name { type Output = Self; + #[inline] fn pow(self, rhs: U256) -> Self::Output { Self::new(self.1.pow(&rhs).retrieve()) } @@ -360,6 +371,7 @@ macro_rules! define_finite_prime_field { /// which we unwrap. Otherwise, there will be panic. impl Rem for $wrapper_name { type Output = Self; + #[inline] fn rem(self, other: Self) -> Self::Output { // create our own check for zeroness? Self::new( @@ -370,6 +382,7 @@ macro_rules! define_finite_prime_field { } } impl Euclid for $wrapper_name { + #[inline] fn div_euclid(&self, other: &Self) -> Self { if other.is_zero() { return Self::from(0u64); @@ -385,6 +398,7 @@ macro_rules! define_finite_prime_field { } Self::new(_q) } + #[inline] fn rem_euclid(&self, other: &Self) -> Self { if other.is_zero() { return Self::from(0u64); @@ -456,6 +470,7 @@ impl Fp { /// function is inherently expensive, and we never call it on the base field, but if /// we did, it's only defined for p=1. Specialized versions exist for all extensions which /// will require the frobenius transformation. + #[inline(always)] pub fn frobenius(&self, exponent: usize) -> Self { match exponent { 1 => self.pow(BN254_FP_MODULUS.value()), @@ -470,6 +485,7 @@ impl Fp { /// prime that is congruent to 3 mod 4. In this case, the sqrt only has the /// possible solution of $\pm pow(n, \frac{p+1}{4})$, which is where this magic /// number below comes from ;) + #[inline] pub fn sqrt(&self) -> CtOption { let arg = ((Self::new(Self::characteristic()) + Self::one()) / Self::from(4)).value(); let sqrt = self.pow(arg); @@ -477,6 +493,7 @@ impl Fp { CtOption::new(sqrt, sqrt.square().ct_eq(self)) } /// Returns the square of the element in the base field + #[inline] pub fn square(&self) -> Self { (*self) * (*self) } @@ -491,7 +508,7 @@ impl Fp { /// Determines the 'sign' of a value in the base field, /// see for more details pub fn sgn0(&self) -> Choice { - let a = *self % Self::from(2u64); + let a = *self % Self::TWO; tracing::debug!(?a, "Fp::sgn0"); if a.is_zero() { Choice::from(0u8) @@ -499,6 +516,26 @@ impl Fp { Choice::from(1u8) } } + /// There is a need to at times move to a representation of the field element with + /// a lower Hamming weight, for instance in the case of multiplication of a group element by + /// such a scalar. This implements the prodinger algorithm, and returns a string of the + /// positive bits and a string of negative bits for the NAF representation + /// see + pub(crate) fn compute_naf(self) -> (U256, U256) { + let x = self.value(); + let xh = x >> 1; + let x3 = x + xh; + let c = xh ^ x3; + let np = x3 & c; + let nm = xh & c; + + (np, nm) + } +} +impl Fr { + pub(crate) fn compute_naf(self) -> (U256, U256) { + Fp::from(self).compute_naf() + } } /// the code below makes the base field "visible" to higher /// order extensions. The issue is really the fact that generic @@ -510,9 +547,6 @@ impl Fp { /// by manually specifying the traits D, N. This enforces the logic /// by means of manual input. impl FieldExtensionTrait<2, 2> for Fp { - fn quadratic_non_residue() -> Self { - >::quadratic_non_residue() - } fn rand(rng: &mut R) -> Self { >::rand(rng) } diff --git a/src/fields/fp12.rs b/src/fields/fp12.rs index dc4b970..a9349ca 100644 --- a/src/fields/fp12.rs +++ b/src/fields/fp12.rs @@ -162,27 +162,11 @@ const FROBENIUS_COEFF_FP12_C1: &[Fp2; 12] = &[ ])), ]), ]; -const FP12_QUADRATIC_NON_RESIDUE: Fp12 = Fp12::new(&[ - Fp6::new(&[ - Fp2::new(&[Fp::ZERO, Fp::ZERO]), - Fp2::new(&[Fp::ZERO, Fp::ZERO]), - Fp2::new(&[Fp::ZERO, Fp::ZERO]), - ]), - Fp6::new(&[ - Fp2::new(&[Fp::ONE, Fp::ZERO]), - Fp2::new(&[Fp::ZERO, Fp::ZERO]), - Fp2::new(&[Fp::ZERO, Fp::ZERO]), - ]), -]); /// type alias for dodectic extension in the representation a + bw for a,b\in Fp6 pub type Fp12 = FieldExtension<12, 2, Fp6>; impl FieldExtensionTrait<12, 2> for Fp12 { - fn quadratic_non_residue() -> Self { - // Self::new(&[Fp6::zero(), Fp6::one()]) - FP12_QUADRATIC_NON_RESIDUE - } fn rand(rng: &mut R) -> Self { Self([ >::rand(rng), @@ -196,10 +180,10 @@ impl FieldExtensionTrait<12, 2> for Fp12 { impl<'a, 'b> Mul<&'b Fp12> for &'a Fp12 { type Output = Fp12; + #[inline] fn mul(self, other: &'b Fp12) -> Self::Output { - // this is again simple Karatsuba multiplication - // see comments in Fp2 impl of `Mul` trait, or otherwise see Alg 20 of - // + // this is simple FOIL'ing of the multiplication of the Fp12 elements in their (Fp6, Fp6) + // representation, see Alg 20 of let t0 = self.0[0] * other.0[0]; let t1 = self.0[1] * other.0[1]; tracing::debug!(?t0, ?t1, "Fp12::mul"); @@ -212,17 +196,20 @@ impl<'a, 'b> Mul<&'b Fp12> for &'a Fp12 { } impl Mul for Fp12 { type Output = Self; + #[inline] fn mul(self, other: Self) -> Self::Output { (&self).mul(&other) } } impl MulAssign for Fp12 { + #[inline] fn mul_assign(&mut self, other: Self) { *self = *self * other; } } impl Inv for Fp12 { type Output = Self; + #[inline] fn inv(self) -> Self::Output { // Implements Alg 23 of let tmp = (self.0[0].square() - (self.0[1].square().residue_mul())).inv(); @@ -232,6 +219,7 @@ impl Inv for Fp12 { } impl One for Fp12 { + #[inline] fn one() -> Self { Self::new(&[Fp6::one(), Fp6::zero()]) } @@ -243,11 +231,13 @@ impl One for Fp12 { #[allow(clippy::suspicious_arithmetic_impl)] impl Div for Fp12 { type Output = Self; + #[inline] fn div(self, other: Self) -> Self::Output { self * other.inv() } } impl DivAssign for Fp12 { + #[inline] fn div_assign(&mut self, other: Self) { *self = *self / other; } @@ -264,6 +254,7 @@ impl ConditionallySelectable for Fp12 { } /// Below are additional functions needed on Fp12 for the pairing operations impl Fp12 { + #[inline] pub(crate) fn unitary_inverse(&self) -> Self { Self::new(&[self.0[0], -self.0[1]]) } @@ -378,6 +369,7 @@ impl Fp12 { Fp12::new(&[Fp6::new(&[z0, z1, z2]), Fp6::new(&[z3, z4, z5])]) } + #[inline(always)] pub(crate) fn frobenius(&self, exponent: usize) -> Self { Self::new(&[ self.0[0].frobenius(exponent), @@ -386,6 +378,7 @@ impl Fp12 { .scale(FROBENIUS_COEFF_FP12_C1[exponent % 12]), ]) } + #[inline] pub(crate) fn square(&self) -> Self { // For F_{p^{12}} = F_{p^6}(w)/(w^2-\gamma), and A=a_0 + a_1*w \in F_{p^{12}}, // we determine C=c_0+c_1*w = A^2\in F_{p^{12}} diff --git a/src/fields/fp2.rs b/src/fields/fp2.rs index 3b1fdce..fc8ad94 100644 --- a/src/fields/fp2.rs +++ b/src/fields/fp2.rs @@ -3,13 +3,12 @@ //! that elements of this field are represented as a_0 + a_1*X. This implements //! the specific behaviour for this extension, such as multiplication. use crate::fields::extensions::FieldExtension; -use crate::fields::fp::{FieldExtensionTrait, Fp, BN254_FP_MODULUS}; +use crate::fields::fp::{FieldExtensionTrait, Fp, BN254_FP_MODULUS, FP_QUADRATIC_NON_RESIDUE}; use crypto_bigint::{rand_core::CryptoRngCore, subtle::ConditionallySelectable, U256}; use num_traits::{Inv, One, Pow, Zero}; use std::ops::{Div, DivAssign, Mul, MulAssign}; use subtle::{Choice, ConstantTimeEq, CtOption}; -const FP2_QUADRATIC_NON_RESIDUE: Fp2 = Fp2::new(&[Fp::NINE, Fp::ONE]); pub(crate) const TWO_INV: Fp = Fp::new(U256::from_words([ 11389680472494603940, 14681934109093717318, @@ -52,7 +51,7 @@ pub type Fp2 = FieldExtension<2, 2, Fp>; // helper functions for us on this specific extension, but // don't generalize to any extension. impl Fp2 { - /// A simple square and multipy algorithm for exponentiation + /// A simple square and multiply algorithm for exponentiation /// # Arguments /// * `by` - Fp, the exponent to raise the element to /// @@ -72,19 +71,26 @@ impl Fp2 { } res } - + #[inline(always)] pub(crate) fn residue_mul(&self) -> Self { - self * &FP2_QUADRATIC_NON_RESIDUE + // Instead of simply `self * &FP2_QUADRATIC_NON_RESIDUE`, we do + // the multiplication "manually", namely: + // (a+bu)*(9+u) = (9a-b)+(a+9b)u, which is cheaper arithmetic in Fp that multiplication + Self::new(&[ + Fp::NINE * self.0[0] - self.0[1], + self.0[0] + Fp::NINE * self.0[1], + ]) } /// Frobenius mapping of a quadratic extension element to a given power /// # Arguments /// * `exponent` - usize, the power to raise the element to + #[inline(always)] pub fn frobenius(&self, exponent: usize) -> Self { let frobenius_coeff_fp2: &[Fp; 2] = &[ // Fp::quadratic_non_residue()**(((p^0) - 1) / 2) Fp::ONE, // Fp::quadratic_non_residue()**(((p^1) - 1) / 2) - >::quadratic_non_residue(), + FP_QUADRATIC_NON_RESIDUE, ]; match exponent % 2 { 0 => *self, @@ -115,15 +121,12 @@ impl Fp2 { } } pub fn square(&self) -> Self { - let t0 = self.0[0] * self.0[1]; - tracing::debug!(?t0, "Fp2::square"); - Self([ - (self.0[1] * >::quadratic_non_residue() + self.0[0]) - * (self.0[0] + self.0[1]) - - t0 - - t0 * >::quadratic_non_residue(), - t0 + t0, - ]) + // We implement manual squaring here and avoid multiplications at all costs + let a = self.0[0] + self.0[1]; + let b = self.0[0] - self.0[1]; + let c = self.0[0] + self.0[0]; + tracing::debug!(?a, "Fp2::square"); + Self([a * b, c * self.0[1]]) } pub fn is_square(&self) -> Choice { let legendre = |x: &Fp| -> i32 { @@ -137,8 +140,7 @@ impl Fp2 { -1 } }; - let sum = self.0[0].square() - + >::quadratic_non_residue() * (-self.0[0]).square(); + let sum = self.0[0].square() + FP_QUADRATIC_NON_RESIDUE * (-self.0[0]).square(); tracing::debug!(?sum, "Fp2::is_square"); Choice::from((legendre(&sum) != -1) as u8) } @@ -150,9 +152,6 @@ impl Fp2 { } } impl FieldExtensionTrait<2, 2> for Fp2 { - fn quadratic_non_residue() -> Self { - FP2_QUADRATIC_NON_RESIDUE - } fn rand(rng: &mut R) -> Self { Self([ >::rand(rng), @@ -167,23 +166,33 @@ impl FieldExtensionTrait<2, 2> for Fp2 { } impl<'a, 'b> Mul<&'b Fp2> for &'a Fp2 { type Output = Fp2; + #[inline] fn mul(self, other: &'b Fp2) -> Self::Output { // This requires a bit more consideration. In Fp2, - // in order to multiply, we must implement complex Karatsuba - // multiplication. - // See https://eprint.iacr.org/2006/471.pdf, Sec 3 - // We create the addition chain from Algo 1 of https://eprint.iacr.org/2022/367.pdf - let t0 = self.0[0] * other.0[0]; - let t1 = self.0[1] * other.0[1]; - + // in order to multiply, we could implement complex Karatsuba + // multiplication, see , Sec 3 + // and then use the addition chain from Alg 1 of : + // // let t0 = self.0[0] * other.0[0]; + // // let t1 = self.0[1] * other.0[1]; + // // + // // Self::Output::new(&[ + // // t1 * FP_QUADRATIC_NON_RESIDUE + t0, + // // (self.0[0] + self.0[1]) * (other.0[0] + other.0[1]) - t0 - t1, + // // ]) + // BUT this is not constant-time, and turns out slower than not invoking the quadratic residue and + // simply doing the schoolbook version, known as the sum of products approach. There is + // an optimized version of this implementation given in Alg 2 of the above reference, + // which requires more granular control of the limb arithmetic over the multiprecision + // scalars than is convenient to implement here. Self::Output::new(&[ - t1 * >::quadratic_non_residue() + t0, - (self.0[0] + self.0[1]) * (other.0[0] + other.0[1]) - t0 - t1, + self.0[0] * other.0[0] - self.0[1] * other.0[1], + self.0[0] * other.0[1] + self.0[1] * other.0[0], ]) } } impl Mul for Fp2 { type Output = Self; + #[inline] fn mul(self, other: Fp2) -> Self::Output { // TODO linter complains about this being a needless reference if I do &a * &b, so this // gets around it @@ -191,6 +200,7 @@ impl Mul for Fp2 { } } impl MulAssign for Fp2 { + #[inline] fn mul_assign(&mut self, other: Self) { *self = *self * other; } @@ -198,12 +208,11 @@ impl MulAssign for Fp2 { impl Inv for Fp2 { type Output = Self; + #[inline] fn inv(self) -> Self { let c0_squared = self.0[0].square(); let c1_squared = self.0[1].square(); - let tmp = (c0_squared - - (c1_squared * >::quadratic_non_residue())) - .inv(); + let tmp = (c0_squared - (c1_squared * FP_QUADRATIC_NON_RESIDUE)).inv(); Self::new(&[self.0[0] * tmp, -(self.0[1] * tmp)]) } } @@ -212,6 +221,7 @@ impl Inv for Fp2 { // this must be defined only for the specific case here, aka not // in extensions.rs impl One for Fp2 { + #[inline] fn one() -> Self { Self::new(&[Fp::ONE, Fp::ZERO]) } @@ -223,11 +233,13 @@ impl One for Fp2 { #[allow(clippy::suspicious_arithmetic_impl)] impl Div for Fp2 { type Output = Self; + #[inline] fn div(self, other: Self) -> Self { self * other.inv() } } impl DivAssign for Fp2 { + #[inline] fn div_assign(&mut self, other: Self) { *self = *self / other; } @@ -246,9 +258,6 @@ impl ConditionallySelectable for Fp2 { // the below is again to make the quadratic extension visible to // higher order sextic extension impl FieldExtensionTrait<6, 3> for Fp2 { - fn quadratic_non_residue() -> Self { - FP2_QUADRATIC_NON_RESIDUE - } fn rand(rng: &mut R) -> Self { >::rand(rng) } @@ -264,6 +273,7 @@ impl FieldExtensionTrait<6, 3> for Fp2 { mod tests { use super::*; use crypto_bigint::U256; + const FP2_QUADRATIC_NON_RESIDUE: Fp2 = Fp2::new(&[Fp::NINE, Fp::ONE]); fn create_field(value: [u64; 4]) -> Fp { Fp::new(U256::from_words(value)) diff --git a/src/fields/fp6.rs b/src/fields/fp6.rs index 9392eab..27a87e6 100644 --- a/src/fields/fp6.rs +++ b/src/fields/fp6.rs @@ -165,18 +165,16 @@ const FROBENIUS_COEFF_FP6_C2: &[Fp2; 6] = &[ ])), ]), ]; -const FP6_QUADRATIC_NON_RESIDUE: Fp6 = Fp6::new(&[ - Fp2::new(&[Fp::ZERO, Fp::ZERO]), - Fp2::new(&[Fp::ONE, Fp::ZERO]), - Fp2::new(&[Fp::ZERO, Fp::ZERO]), -]); + /// type alias for the sextic extension of the base field pub type Fp6 = FieldExtension<6, 3, Fp2>; impl Fp6 { + #[inline(always)] pub(crate) fn residue_mul(&self) -> Self { Self([self.0[2].residue_mul(), self.0[0], self.0[1]]) } + #[inline(always)] pub(crate) fn frobenius(&self, exponent: usize) -> Self { Self::new(&[ self.0[0].frobenius(exponent), @@ -209,10 +207,6 @@ impl Fp6 { } } impl FieldExtensionTrait<6, 3> for Fp6 { - fn quadratic_non_residue() -> Self { - FP6_QUADRATIC_NON_RESIDUE - } - fn rand(rng: &mut R) -> Self { Self([ >::rand(rng), @@ -226,28 +220,102 @@ impl FieldExtensionTrait<6, 3> for Fp6 { } impl<'a, 'b> Mul<&'b Fp6> for &'a Fp6 { type Output = Fp6; + #[inline] fn mul(self, other: &'b Fp6) -> Self::Output { - // This is the exact same strategy as multiplication in Fp2 - // see the doc string there for more details - let t0 = self.0[0] * other.0[0]; - let t1 = self.0[1] * other.0[1]; - let t2 = self.0[2] * other.0[2]; - tracing::debug!(?t0, ?t1, ?t2, "Fp6::mul"); + // We could do Karatsuba multiplication here, which would look simpler: + // // let t0 = self.0[0] * other.0[0]; + // // let t1 = self.0[1] * other.0[1]; + // // let t2 = self.0[2] * other.0[2]; + // // tracing::debug!(?t0, ?t1, ?t2, "Fp6::mul"); + // // + // // Self::Output::new(&[ + // // ((self.0[1] + self.0[2]) * (other.0[1] + other.0[2]) - t1 - t2).residue_mul() + t0, + // // (self.0[0] + self.0[1]) * (other.0[0] + other.0[1]) - t0 - t1 + t2.residue_mul(), + // // (self.0[0] + self.0[2]) * (other.0[0] + other.0[2]) - t0 + t1 - t2, + // // ]) + // + // But the issue again is constant-time execution. We opt for schoolbook multiplication + // here instead following Algo 5 of , which yields + // the following results: + // + // c0,0 = a0,0b0,0 - a0,1b0,1 + αa1,0b2,0 - αa1,1b2,1 + αa2,0b1,0 - αa2,1b1,1 - a1,0b2,1 - a1,1b2,0 + // - a2,0b1,1 - a2,1b1,0. + // = a0,0b0,0 - a0,1b0,1 + a1,0(αb2,0 - b2,1) - a1,1(b2,0 + αb2,1) + a2,0(αb1,0 - b1,1) + // - a2,1(b1,0 + αb1,1). + // c0,1 = a0,0b0,1 + a0,1b0,0 + αa1,0b2,1 + αa1,1b2,0 + αa2,0b1,1 + αa2,1b1,0 + a1,0b2,0 - a1,1b2,1 + // + a2,0b1,0 - a2,1b1,1. + // = a0,0b0,1 + a0,1b0,0 + a1,0(b2,0 + αb2,1) + a1,1(αb2,0 - b2,1) + a2,0(b1,0 + αb1,1) + // + a2,1(αb1,0 - b1,1). + // c1,0 = a0,0b1,0 - a0,1b1,1 + a1,0b0,0 - a1,1b0,1 + αa2,0b2,0 - αa2,1b2,1 - a2,0b2,1 - a2,1b2,0. + // = a0,0b1,0 - a0,1b1,1 + a1,0b0,0 - a1,1b0,1 + a2,0(αb2,0 - b2,1) - a2,1(αb2,1 + b2,0). + // c1,1 = a0,0b1,1 + a0,1b1,0 + a1,0b0,1 + a1,1b0,0 + αa2,0b2,1 + αa2,1b2,0 + a2,0b2,0 - a2,1b2,1. + // = a0,0b1,1 + a0,1b1,0 + a1,0b0,1 + a1,1b0,0 + a2,0(αb2,1 + b2,0) + a2,1(αb2,0 - b2,1). + // c2,0 = a0,0b2,0 - a0,1b2,1 + a1,0b1,0 - a1,1b1,1 + a2,0b0,0 - a2,1b0,1. + // c2,1 = a0,0b2,1 + a0,1b2,0 + a1,0b1,1 + a1,1b1,0 + a2,0b0,1 + a2,1b0,0. + // + // Here, alpha is \xi=\alpha+u=9+u => alpha = 9 + + let a20_m_b21 = Fp::NINE * other.0[2].0[0] - other.0[2].0[1]; + let a10_m_b11 = Fp::NINE * other.0[1].0[0] - other.0[1].0[1]; + let b21_p_b20 = Fp::NINE * other.0[2].0[1] + other.0[2].0[0]; + let b20_m_b21 = Fp::NINE * other.0[2].0[0] - other.0[2].0[1]; + let b11_p_b10 = Fp::NINE * other.0[1].0[1] + other.0[1].0[0]; + + let c00 = self.0[0].0[0] * other.0[0].0[0] - self.0[0].0[1] * other.0[0].0[1] + + self.0[1].0[0] * a20_m_b21 + - self.0[1].0[1] * b21_p_b20 + + self.0[2].0[0] * a10_m_b11 + - self.0[2].0[1] * b11_p_b10; + + let c01 = self.0[0].0[0] * other.0[0].0[1] + + self.0[0].0[1] * other.0[0].0[0] + + self.0[1].0[0] * b21_p_b20 + + self.0[1].0[1] * b20_m_b21 + + self.0[2].0[0] * b11_p_b10 + + self.0[2].0[1] * a10_m_b11; + + let c10 = self.0[0].0[0] * other.0[1].0[0] - self.0[0].0[1] * other.0[1].0[1] + + self.0[1].0[0] * other.0[0].0[0] + - self.0[1].0[1] * other.0[0].0[1] + + self.0[2].0[0] * b20_m_b21 + - self.0[2].0[1] * b21_p_b20; + + let c11 = self.0[0].0[0] * other.0[1].0[1] + + self.0[0].0[1] * other.0[1].0[0] + + self.0[1].0[0] * other.0[0].0[1] + + self.0[1].0[1] * other.0[0].0[0] + + self.0[2].0[0] * b21_p_b20 + + self.0[2].0[1] * a20_m_b21; + + let c20 = self.0[0].0[0] * other.0[2].0[0] - self.0[0].0[1] * other.0[2].0[1] + + self.0[1].0[0] * other.0[1].0[0] + - self.0[1].0[1] * other.0[1].0[1] + + self.0[2].0[0] * other.0[0].0[0] + - self.0[2].0[1] * other.0[0].0[1]; + + let c21 = self.0[0].0[0] * other.0[2].0[1] + + self.0[0].0[1] * other.0[2].0[0] + + self.0[1].0[0] * other.0[1].0[1] + + self.0[1].0[1] * other.0[1].0[0] + + self.0[2].0[0] * other.0[0].0[1] + + self.0[2].0[1] * other.0[0].0[0]; Self::Output::new(&[ - ((self.0[1] + self.0[2]) * (other.0[1] + other.0[2]) - t1 - t2).residue_mul() + t0, - (self.0[0] + self.0[1]) * (other.0[0] + other.0[1]) - t0 - t1 + t2.residue_mul(), - (self.0[0] + self.0[2]) * (other.0[0] + other.0[2]) - t0 + t1 - t2, + Fp2::new(&[c00, c01]), + Fp2::new(&[c10, c11]), + Fp2::new(&[c20, c21]), ]) } } impl Mul for Fp6 { type Output = Self; + #[inline] fn mul(self, other: Self) -> Self::Output { (&self).mul(&other) } } impl MulAssign for Fp6 { + #[inline] fn mul_assign(&mut self, other: Self) { *self = *self * other; } @@ -255,6 +323,7 @@ impl MulAssign for Fp6 { impl Inv for Fp6 { type Output = Self; + #[inline] fn inv(self) -> Self::Output { // Implements a low-overhead version of Alg 17 of let t0 = self.0[0].square() - self.0[1] * self.0[2].residue_mul(); @@ -268,6 +337,7 @@ impl Inv for Fp6 { } impl One for Fp6 { + #[inline] fn one() -> Self { Self::new(&[Fp2::one(), Fp2::zero(), Fp2::zero()]) } @@ -279,11 +349,13 @@ impl One for Fp6 { #[allow(clippy::suspicious_arithmetic_impl)] impl Div for Fp6 { type Output = Self; + #[inline] fn div(self, other: Self) -> Self::Output { self * other.inv() } } impl DivAssign for Fp6 { + #[inline] fn div_assign(&mut self, other: Self) { *self = *self / other; } @@ -302,9 +374,6 @@ impl ConditionallySelectable for Fp6 { // make sextic extension visible to the dodectic extension impl FieldExtensionTrait<12, 2> for Fp6 { - fn quadratic_non_residue() -> Self { - >::quadratic_non_residue() - } fn rand(rng: &mut R) -> Self { >::rand(rng) } diff --git a/src/groups/g1.rs b/src/groups/g1.rs index 8a46359..8b0f8c7 100644 --- a/src/groups/g1.rs +++ b/src/groups/g1.rs @@ -82,35 +82,24 @@ impl G1Affine { // # Arguments // * `v` - a tuple of field elements that represent the x and y coordinates of the point fn new(v: [Fp; 2]) -> Result { - let _g1affine_is_on_curve = |x: &Fp, y: &Fp, z: &Choice| -> Choice { - let y2 = y.square(); - let x2 = x.square(); - let lhs = y2 - (x2 * (*x)); + let is_on_curve = { + let y2 = v[1].square(); + let x2 = v[0].square(); + let lhs = y2 - (x2 * v[0]); let rhs = >::curve_constant(); tracing::debug!(?y2, ?x2, ?lhs, ?rhs, "G1Affine::new"); - lhs.ct_eq(&rhs) | *z + lhs.ct_eq(&rhs) }; - let _g1affine_is_torsion_free = |_x: &Fp, _y: &Fp, _z: &Choice| -> Choice { - // every point in G1 on the curve is in the r-torsion of BN254 - Choice::from(1u8) - }; - let is_on_curve: Choice = _g1affine_is_on_curve(&v[0], &v[1], &Choice::from(0u8)); + // every point in G1 on the curve is in the r-torsion of BN254, + // so we don't need to check for subgroup membership tracing::debug!(?is_on_curve, "G1Affine::new"); match bool::from(is_on_curve) { - true => { - let is_in_torsion: Choice = - _g1affine_is_torsion_free(&v[0], &v[1], &Choice::from(0u8)); - tracing::debug!(?is_in_torsion, "G1Affine::new"); - match bool::from(is_in_torsion) { - true => Ok(Self { - x: v[0], - y: v[1], - infinity: Choice::from(0u8), - }), - _ => Err(GroupError::NotInSubgroup), - } - } + true => Ok(Self { + x: v[0], + y: v[1], + infinity: Choice::from(0u8), + }), false => Err(GroupError::NotOnCurve), } } @@ -181,32 +170,22 @@ impl G1Projective { /// * `v` - a tuple of field elements that represent the x, y, and z coordinates of the point #[allow(dead_code)] pub fn new(v: [Fp; 3]) -> Result { - let _g1projective_is_on_curve = |x: &Fp, y: &Fp, z: &Fp| -> Choice { - let y2 = y.square(); - let x2 = x.square(); - let z2 = z.square(); - let lhs = y2 * (*z); - let rhs = x2 * (*x) + z2 * (*z) * >::curve_constant(); + let is_on_curve = { + let y2 = v[1].square(); + let x2 = v[0].square(); + let z2 = v[2].square(); + let lhs = y2 * v[2]; + let rhs = x2 * v[0] + z2 * v[2] * >::curve_constant(); tracing::debug!(?y2, ?x2, ?z2, ?lhs, ?rhs, "G1Projective::new"); - lhs.ct_eq(&rhs) | Choice::from(z.is_zero() as u8) + lhs.ct_eq(&rhs) | Choice::from(v[2].is_zero() as u8) }; - let _g1projective_is_torsion_free = - |_x: &Fp, _y: &Fp, _z: &Fp| -> Choice { Choice::from(1u8) }; - let is_on_curve: Choice = _g1projective_is_on_curve(&v[0], &v[1], &v[2]); tracing::debug!(?is_on_curve, "G1Projective::new"); match bool::from(is_on_curve) { - true => { - let is_in_torsion: Choice = _g1projective_is_torsion_free(&v[0], &v[1], &v[2]); - tracing::debug!(?is_in_torsion, "G1Projective::new"); - match bool::from(is_in_torsion) { - true => Ok(Self { - x: v[0], - y: v[1], - z: v[2], - }), - false => Err(GroupError::NotOnCurve), - } - } + true => Ok(Self { + x: v[0], + y: v[1], + z: v[2], + }), false => Err(GroupError::NotOnCurve), } } diff --git a/src/groups/g2.rs b/src/groups/g2.rs index 45f22a3..e74ca1b 100644 --- a/src/groups/g2.rs +++ b/src/groups/g2.rs @@ -219,16 +219,15 @@ impl G2Affine { /// # Arguments /// * `v` - a tuple of field elements that represent the x and y coordinates of the point fn new_unchecked(v: [Fp2; 2]) -> Result { - let _g2affine_is_on_curve = |x: &Fp2, y: &Fp2, z: &Choice| -> Choice { - let y2 = y.square(); - let x2 = x.square(); - let lhs = y2 - (x2 * (*x)); + let is_on_curve = { + let y2 = v[1].square(); + let x2 = v[0].square(); + let lhs = y2 - (x2 * v[0]); let rhs = >::curve_constant(); tracing::debug!(?y2, ?x2, ?lhs, ?rhs, "G2Affine::new_unchecked"); - lhs.ct_eq(&rhs) | *z + lhs.ct_eq(&rhs) }; - let is_on_curve = _g2affine_is_on_curve(&v[0], &v[1], &Choice::from(0u8)); match bool::from(is_on_curve) { true => Ok(Self { x: v[0], @@ -247,14 +246,14 @@ impl G2Projective { /// # Arguments /// * `v` - a tuple of field elements that represent the x, y, and z coordinates of the point pub fn new(v: [Fp2; 3]) -> Result { - let _g2projective_is_on_curve = |x: &Fp2, y: &Fp2, z: &Fp2| -> Choice { - let y2 = y.square(); - let x2 = x.square(); - let z2 = z.square(); - let lhs = y2 * (*z); - let rhs = x2 * (*x) + z2 * (*z) * >::curve_constant(); + let is_on_curve = { + let y2 = v[1].square(); + let x2 = v[0].square(); + let z2 = v[2].square(); + let lhs = y2 * v[2]; + let rhs = x2 * v[0] + z2 * v[2] * >::curve_constant(); tracing::debug!(?y2, ?x2, ?z2, ?lhs, ?rhs, "G2Projective::new"); - lhs.ct_eq(&rhs) | Choice::from(z.is_zero() as u8) + lhs.ct_eq(&rhs) | Choice::from(v[2].is_zero() as u8) }; // This method is where the magic happens. In a naïve approach, in order to check for // validity in the r-torsion, one could simply verify the r-torsion condition: @@ -279,11 +278,11 @@ impl G2Projective { // ---------- // 1. // 2. - let _g2projective_is_torsion_free = |x: &Fp2, y: &Fp2, z: &Fp2| -> Choice { + let is_torsion_free = { let tmp = G2Projective { - x: *x, - y: *y, - z: *z, + x: v[0], + y: v[1], + z: v[2], }; let mut a = tmp * BLS_X; // xQ let b = a.endomorphism(); // ψ(xQ) @@ -292,9 +291,7 @@ impl G2Projective { let lhs = rhs + b + a; // ψ^2(xQ) + ψ(xQ) + (x+1)Q rhs = rhs.endomorphism().double() - lhs; // ψ^3(2xQ) - (ψ^2(xQ) + ψ(xQ) + (x+1)Q) tracing::debug!( - ?x, - ?y, - ?z, + ?v, ?a, ?b, ?lhs, @@ -305,10 +302,8 @@ impl G2Projective { // we do two checks: one is to verify that the result is indeed a point at infinity, // but we need a second check to verify that it is OUR point at infinity, namely for // the curve defined on the twist. - Choice::from(rhs.is_zero() as u8) & _g2projective_is_on_curve(&rhs.x, &rhs.y, &rhs.z) + Choice::from(rhs.is_zero() as u8) & is_on_curve }; - let is_on_curve = _g2projective_is_on_curve(&v[0], &v[1], &v[2]); - let is_torsion_free = _g2projective_is_torsion_free(&v[0], &v[1], &v[2]); match bool::from(is_on_curve) { true => match bool::from(is_torsion_free) { true => Ok(Self { diff --git a/src/groups/group.rs b/src/groups/group.rs index ec04fcd..11a905f 100644 --- a/src/groups/group.rs +++ b/src/groups/group.rs @@ -146,7 +146,7 @@ impl<'a, const D: usize, const N: usize, F: FieldExtensionTrait> Neg for &'a GroupAffine { type Output = GroupAffine; - + #[inline] fn neg(self) -> Self::Output { Self::Output { x: self.x, @@ -157,6 +157,7 @@ impl<'a, const D: usize, const N: usize, F: FieldExtensionTrait> Neg } impl> Neg for GroupAffine { type Output = GroupAffine; + #[inline] fn neg(self) -> Self::Output { -&self } @@ -190,6 +191,7 @@ impl> Conditionally impl> PartialEq for GroupAffine { + #[inline] fn eq(&self, other: &Self) -> bool { bool::from(self.ct_eq(other)) } @@ -202,7 +204,7 @@ impl> GroupAffine bool { bool::from(self.infinity) } @@ -235,7 +237,7 @@ impl> GroupProjecti z: F::zero(), } } - + #[inline(always)] pub(crate) fn is_zero(&self) -> bool { self.z.is_zero() } @@ -294,7 +296,7 @@ impl<'a, const D: usize, const N: usize, F: FieldExtensionTrait> Neg for &'a GroupProjective { type Output = GroupProjective; - + #[inline] fn neg(self) -> Self::Output { Self::Output { x: self.x, @@ -307,6 +309,7 @@ impl> Neg for GroupProjective { type Output = Self; + #[inline] fn neg(self) -> Self::Output { -&self } @@ -347,6 +350,7 @@ impl> Conditionally impl> PartialEq for GroupProjective { + #[inline] fn eq(&self, other: &Self) -> bool { bool::from(self.ct_eq(other)) } @@ -411,6 +415,7 @@ impl<'a, 'b, const D: usize, const N: usize, F: FieldExtensionTrait> /// /// Complexity: /// `12M` + `2m` + `19A` + #[inline] fn add(self, other: &'b GroupProjective) -> Self::Output { let t0 = self.x * other.x; let t1 = self.y * other.y; @@ -476,6 +481,7 @@ impl> Add { type Output = Self; + #[inline] fn add(self, rhs: GroupProjective) -> Self::Output { &self + &rhs } @@ -486,6 +492,7 @@ impl<'a, 'b, const D: usize, const N: usize, F: FieldExtensionTrait> Sub<&'b GroupProjective> for &'a GroupProjective { type Output = GroupProjective; + #[inline] fn sub(self, other: &'b GroupProjective) -> Self::Output { self + &(-other) } @@ -494,6 +501,7 @@ impl> Sub { type Output = Self; + #[inline] fn sub(self, rhs: GroupProjective) -> Self::Output { &self - &rhs } @@ -504,19 +512,25 @@ impl<'a, 'b, const D: usize, const N: usize, F: FieldExtensionTrait> Mul<& for &'a GroupProjective { /// This is simply the `double-and-add` algorithm for multiplication, which is the ECC - /// equivalent of the `square-and-multiply` algorithm used in modular exponentiation. + /// equivalent of the `square-and-multiply` algorithm used in modular exponentiation. It uses + /// the lower Hamming weight representation of the scalar to reduce the number of operations /// /// type Output = GroupProjective; fn mul(self, other: &'b Fp) -> Self::Output { - let bits = other.value().to_le_bytes(); + let (np, nm) = other.compute_naf(); let mut res = Self::Output::zero(); - for bit in bits.iter().rev() { - for i in (0..8).rev() { - res = res.double(); - if (bit & (1 << i)) != 0 { - res = &res + self; - } + + for i in (0..256).rev() { + res = res.double(); + + let np_bit = np.bit(i).into(); + let nm_bit = nm.bit(i).into(); + + if np_bit { + res = &res + self; + } else if nm_bit { + res = &res - self; } } res @@ -527,6 +541,7 @@ impl> Mul for GroupProjective { type Output = Self; + #[inline] fn mul(self, rhs: Fp) -> Self::Output { &self * &rhs } diff --git a/src/groups/gt.rs b/src/groups/gt.rs index 15c0003..f79578c 100644 --- a/src/groups/gt.rs +++ b/src/groups/gt.rs @@ -140,6 +140,7 @@ impl ConditionallySelectable for Gt { } impl PartialEq for Gt { + #[inline] fn eq(&self, other: &Self) -> bool { bool::from(self.ct_eq(other)) } @@ -166,19 +167,25 @@ impl<'a, 'b> Sub<&'b Gt> for &'a Gt { #[allow(clippy::suspicious_arithmetic_impl)] impl<'a, 'b> Mul<&'b Fr> for &'a Gt { /// This is simply the `double-and-add` algorithm for multiplication, which is the ECC - /// equivalent of the `square-and-multiply` algorithm used in modular exponentiation. + /// equivalent of the `square-and-multiply` algorithm used in modular exponentiation. It uses + // the lower Hamming weight representation of the scalar to reduce the number of operations /// /// type Output = Gt; fn mul(self, other: &'b Fr) -> Self::Output { - let bits = other.value().to_le_bytes(); + let (np, nm) = other.compute_naf(); let mut res = Self::Output::identity(); - for bit in bits.iter().rev() { - for i in (0..8).rev() { - res = res.double(); - if (bit & (1 << i)) != 0 { - res = &res + self; - } + + for i in (0..256).rev() { + res = res.double(); + + let np_bit = np.bit(i).into(); + let nm_bit = nm.bit(i).into(); + + if np_bit { + res = &res + self; + } else if nm_bit { + res = &res - self; } } res @@ -187,6 +194,7 @@ impl<'a, 'b> Mul<&'b Fr> for &'a Gt { impl Mul for Gt { type Output = Self; + #[inline] fn mul(self, rhs: Fr) -> Self::Output { &self * &rhs } diff --git a/src/lib.rs b/src/lib.rs index 1646cc2..4bd8720 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -2,7 +2,11 @@ #![doc = include_str!("../README.md")] #![deny(unsafe_code)] #![deny(dead_code)] -#![allow(clippy::needless_doctest_main, clippy::doc_lazy_continuation)] +#![allow( + clippy::needless_doctest_main, + clippy::doc_lazy_continuation, + clippy::too_long_first_doc_paragraph +)] #![warn( clippy::unwrap_used, missing_docs, diff --git a/src/pairing.rs b/src/pairing.rs index f80a30c..fc34167 100644 --- a/src/pairing.rs +++ b/src/pairing.rs @@ -102,6 +102,7 @@ impl MillerLoopResult { /// This implements algorithm 9 from https://eprint.iacr.org/2010/354.pdf, with the notable /// difference that instead of passing an element of Fp4 (which I did not implement), we pass /// in only the two components from Fp2 that comprise the Fp4 element. + #[must_use] fn fp4_square(a: Fp2, b: Fp2) -> (Fp2, Fp2) { // Line 1 let t0 = a.square(); @@ -121,6 +122,7 @@ impl MillerLoopResult { /// This implements efficient squaring of an element of Fp12 in the cyclotomic subgroup /// C_{\phi^6}. It is what's called "Granger-Scott" squaring, and is an implementation of /// algorithm 5.5.4 (listing 21) from https://www.math.u-bordeaux.fr/~damienrobert/csi/book/book.pdf + #[must_use] fn cyclotomic_squared(f: Fp12) -> Fp12 { // Lines 3-8 let mut z0 = f.0[0].0[0]; @@ -162,6 +164,7 @@ impl MillerLoopResult { /// This is a simple square and multiply algorithm for exponentiation. You can get more /// complicated algorithms if you go to a compressed representation, such as Algorithm /// 5.5.4, listing 27 + #[must_use] pub(crate) fn cyclotomic_exp(f: Fp12, exponent: &Fp) -> Fp12 { let bits = exponent.value().to_words(); let mut res = Fp12::one(); @@ -380,11 +383,7 @@ impl G2Projective { self.y = e * (i - j) - h * self.y; self.z *= h; - Ell( - >::quadratic_non_residue() * (e * base.x - d * base.y), - d, - e.neg(), - ) + Ell((e * base.x - d * base.y).residue_mul(), d, e.neg()) } fn doubling_step(&mut self) -> Ell { let a = (self.x * self.y).scale(TWO_INV); @@ -404,11 +403,7 @@ impl G2Projective { self.y = g.square() - (e_sq + e_sq + e_sq); self.z = b * h; - Ell( - >::quadratic_non_residue() * i, - h.neg(), - j + j + j, - ) + Ell(i.residue_mul(), h.neg(), j + j + j) } } /// Execute the optimal ate pairing on BN254 for a given input pair of (G1,G2) points. If either