ed448-goldilocks: account for oddness in Scalar divisions

ed448-goldilocks/Cargo.toml

@@ -39,6 +39,7 @@ serde = ["dep:serdect", "ed448?/serde_bytes"]
 [dev-dependencies]
 hex-literal = "1"
 hex = "0.4"
+proptest = "1"
 rand_core = { version = "0.9", features = ["os_rng"] }
 rand_chacha = "0.9"
 serde_bare = "0.5"

ed448-goldilocks/src/edwards/extended.rs

@@ -5,7 +5,7 @@ use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
 use crate::curve::scalar_mul::variable_base;
 use crate::curve::twedwards::extended::ExtendedPoint as TwistedExtendedPoint;
-use crate::field::FieldElement;
+use crate::field::{ConstMontyType, FieldElement};
 use crate::*;
 use elliptic_curve::{
     CurveGroup, Error,
@@ -548,34 +548,7 @@ impl EdwardsPoint {
         scalar_div_four.div_by_four();
 
         // Use isogeny and dual isogeny to compute phi^-1((s/4) * phi(P))
-        let partial_result = variable_base(&self.to_twisted(), &scalar_div_four).to_untwisted();
-        // Add partial result to (scalar mod 4) * P
-        partial_result.add(&self.scalar_mod_four(scalar))
-    }
-
-    /// Returns (scalar mod 4) * P in constant time
-    pub(crate) fn scalar_mod_four(&self, scalar: &EdwardsScalar) -> Self {
-        // Compute compute (scalar mod 4)
-        let s_mod_four = scalar[0] & 3;
-
-        // Compute all possible values of (scalar mod 4) * P
-        let zero_p = EdwardsPoint::IDENTITY;
-        let one_p = self;
-        let two_p = one_p.double();
-        let three_p = two_p.add(self);
-
-        // Under the reasonable assumption that `==` is constant time
-        // Then the whole function is constant time.
-        // This should be cheaper than calling double_and_add or a scalar mul operation
-        // as the number of possibilities are so small.
-        // XXX: This claim has not been tested (although it sounds intuitive to me)
-        let mut result = EdwardsPoint::IDENTITY;
-        result.conditional_assign(&zero_p, Choice::from((s_mod_four == 0) as u8));
-        result.conditional_assign(one_p, Choice::from((s_mod_four == 1) as u8));
-        result.conditional_assign(&two_p, Choice::from((s_mod_four == 2) as u8));
-        result.conditional_assign(&three_p, Choice::from((s_mod_four == 3) as u8));
-
-        result
+        variable_base(&self.to_twisted(), &scalar_div_four).to_untwisted()
     }
 
     /// Standard compression; store Y and sign of X
@@ -722,8 +695,57 @@ impl EdwardsPoint {
     ///   prime-order subgroup;
     /// * `false` if `self` has a nonzero torsion component and is not
     ///   in the prime-order subgroup.
+    // See https://eprint.iacr.org/2022/1164.
     pub fn is_torsion_free(&self) -> Choice {
-        (self * EdwardsScalar::new(ORDER)).ct_eq(&Self::IDENTITY)
+        const A: FieldElement = FieldElement(ConstMontyType::new(&U448::from_be_hex(
+            "fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffffffffffffffffffffffffffffffffffffffffffffffffeceaf",
+        )));
+        const A1: FieldElement = FieldElement(ConstMontyType::new(&U448::from_u64(156320)));
+        const MINUS_SQRT_B1: FieldElement = FieldElement(ConstMontyType::new(&U448::from_be_hex(
+            "749a7410536c225f1025ca374176557d7839611d691caad26d74a1fca5cfad15f196642c0a4484b67f321025577cc6b5a6f443c2eaa36327",
+        )));
+
+        let mut e = self.X * (self.Z - self.Y);
+        let ee = e.square();
+        let mut u = FieldElement::A_PLUS_TWO_OVER_FOUR * (self.Z + self.Y) * e * self.X;
+        let w = self.Z.double() * (self.Z - self.Y);
+
+        let u2 = u.double().double();
+        let w2 = w.double();
+
+        let mut w1 = u2.sqrt();
+        let mut ok = w1.square().ct_eq(&u2);
+        let u1 = (u2 - A1 * ee - w1 * w2).half();
+
+        // If `u1` happens not to be a square, then `sqrt(u1)` returns `sqrt(-u1)`
+        // in that case (since we are in a finite field GF(q) with q = 3 mod 4,
+        // if `u1` is not a square then `-u1` must be a square). In such a case, we
+        // should replace `(u1,w1)` with `((B1*e^4)/u1, -w1)`. To avoid the division,
+        // we instead switch to an isomorphic curve; namely:
+        //   u2 = B1*(e^4)*u1
+        //   w2 = -w1*u1
+        //   e2 = e*u1
+        // Then:
+        //   w = sqrt(u2) = sqrt(-B1)*(e^2)*sqrt(-u1)
+        //   u = (w^2 - A*e^2 - w*w1)/2
+        let mut w = u1.sqrt();
+        let u1_is_square = w.square().ct_eq(&u1);
+        w1.conditional_assign(&-(w1 * u1), !u1_is_square);
+        e.conditional_assign(&(e * u1), !u1_is_square);
+        w.conditional_assign(&(MINUS_SQRT_B1 * ee * w), !u1_is_square);
+        u = (w.square() - A * e.square() - w * w1).half();
+
+        ok &= u.is_square();
+
+        // If the source point was a low-order point, then the computations
+        // above are incorrect. We handle this case here; among the
+        // low-order points, only the neutral point is in the prime-order
+        // subgroup.
+        let is_low_order = self.X.is_zero() | self.Y.is_zero();
+        let is_neutral = self.Y.ct_eq(&self.Z);
+        ok ^= is_low_order & (ok ^ is_neutral);
+
+        ok
     }
 
     /// Hash a message to a point on the curve
@@ -972,6 +994,8 @@ mod tests {
     use super::*;
     use elliptic_curve::Field;
     use hex_literal::hex;
+    use proptest::prelude::any;
+    use proptest::proptest;
     use rand_core::TryRngCore;
 
     fn hex_to_field(hex: &'static str) -> FieldElement {
@@ -1267,4 +1291,15 @@ mod tests {
 
         assert_eq!(computed_commitment, expected_commitment);
     }
+
+    proptest! {
+        #[test]
+        fn fuzz_is_torsion_free(
+           bytes in any::<[u8; 57]>()
+        ) {
+            let scalar = EdwardsScalar::from_bytes_mod_order(&bytes.into());
+            let point = EdwardsPoint::mul_by_generator(&scalar);
+            assert_eq!(point.is_torsion_free().unwrap_u8(), 1);
+        }
+    }
 }

ed448-goldilocks/src/field/element.rs

@@ -9,7 +9,7 @@ use crate::{
 use elliptic_curve::{
     array::Array,
     bigint::{
-        Integer, NonZero, U448, U704,
+        Integer, NonZero, U448, U704, Zero,
         consts::{U56, U84, U88},
     },
     group::cofactor::CofactorGroup,
@@ -261,6 +261,10 @@ impl FieldElement {
         self.0.retrieve().is_odd()
     }
 
+    pub fn is_zero(&self) -> Choice {
+        self.0.is_zero()
+    }
+
     /// Inverts a field element
     /// Previous chain length: 462, new length 460
     pub fn invert(&self) -> Self {
@@ -368,6 +372,10 @@ impl FieldElement {
         (inv_sqrt_x * u, zero_u | is_res)
     }
 
+    pub(crate) fn half(&self) -> FieldElement {
+        Self(self.0.div_by_2())
+    }
+
     pub(crate) fn map_to_curve_elligator2(&self) -> AffinePoint {
         let mut t1 = self.square(); // 1.   t1 = u^2
         t1 *= Self::Z; // 2.   t1 = Z * t1              // Z * u^2

ed448-goldilocks/src/field/scalar.rs

@@ -13,7 +13,7 @@ use elliptic_curve::{
         Array, ArraySize,
         typenum::{Prod, Unsigned},
     },
-    bigint::{Limb, NonZero, U448, U896, Word, Zero},
+    bigint::{Integer, Limb, NonZero, U448, U896, Word, Zero},
     consts::U2,
     ff::{Field, helpers},
     ops::{Invert, Reduce, ReduceNonZero},
@@ -670,6 +670,18 @@ impl<C: CurveWithScalar> Scalar<C> {
     /// This is used in the 2-isogeny when mapping points from Ed448-Goldilocks
     /// to Twisted-Goldilocks
     pub(crate) fn div_by_four(&mut self) {
+        let s_mod_4 = self[0] & 3;
+
+        let s_plus_l = self.scalar + ORDER;
+        let s_plus_2l = s_plus_l + ORDER;
+        let s_plus_3l = s_plus_2l + ORDER;
+
+        self.scalar.conditional_assign(&s_plus_l, s_mod_4.ct_eq(&1));
+        self.scalar
+            .conditional_assign(&s_plus_2l, s_mod_4.ct_eq(&2));
+        self.scalar
+            .conditional_assign(&s_plus_3l, s_mod_4.ct_eq(&3));
+
         self.scalar >>= 2;
     }
 
@@ -786,8 +798,12 @@ impl<C: CurveWithScalar> Scalar<C> {
     }
 
     /// Halves a Scalar modulo the prime
-    pub const fn halve(&self) -> Self {
-        Self::new(self.scalar.shr_vartime(1))
+    pub fn halve(&self) -> Self {
+        let is_odd = self.scalar.is_odd();
+        let if_odd = self.scalar + ORDER;
+        let scalar = U448::conditional_select(&self.scalar, &if_odd, is_odd);
+
+        Self::new(scalar >> 1)
     }
 
     /// Attempt to construct a `Scalar` from a canonical byte representation.