curve: Add arbitrary integer multiplication with `MontgomeryPoint::mu…

…l_bits_be` (#555)

There is occasionally [a need](#519 (comment)) to multiply a non-prime-order Montgomery point by an integer. There's currently no way to do this, since our only methods are multiplication by `Scalar` (doesn't make sense in the non-prime-order case), and `MontgomeryPoint::mul_base_clamped` clamps the integer before multiplying.

This defines `MontgomeryPoint::mul_bits_be`, which takes a big-endian representation of an integer and multiplies the point by that integer. Its usage is not recommended by default, but it is also not so unsafe as to be gated behind a `hazmat` feature.

curve25519-dalek/src/montgomery.rs

@@ -151,6 +151,43 @@ impl MontgomeryPoint {
         Self::mul_base(&s)
     }
 
+    /// Given `self` \\( = u\_0(P) \\), and a big-endian bit representation of an integer
+    /// \\(n\\), return \\( u\_0(\[n\]P) \\). This is constant time in the length of `bits`.
+    ///
+    /// **NOTE:** You probably do not want to use this function. Almost every protocol built on
+    /// Curve25519 uses _clamped multiplication_, explained
+    /// [here](https://neilmadden.blog/2020/05/28/whats-the-curve25519-clamping-all-about/).
+    /// When in doubt, use [`Self::mul_clamped`].
+    pub fn mul_bits_be(&self, bits: impl Iterator<Item = bool>) -> MontgomeryPoint {
+        // Algorithm 8 of Costello-Smith 2017
+        let affine_u = FieldElement::from_bytes(&self.0);
+        let mut x0 = ProjectivePoint::identity();
+        let mut x1 = ProjectivePoint {
+            U: affine_u,
+            W: FieldElement::ONE,
+        };
+
+        // Go through the bits from most to least significant, using a sliding window of 2
+        let mut prev_bit = false;
+        for cur_bit in bits {
+            let choice: u8 = (prev_bit ^ cur_bit) as u8;
+
+            debug_assert!(choice == 0 || choice == 1);
+
+            ProjectivePoint::conditional_swap(&mut x0, &mut x1, choice.into());
+            differential_add_and_double(&mut x0, &mut x1, &affine_u);
+
+            prev_bit = cur_bit;
+        }
+        // The final value of prev_bit above is scalar.bits()[0], i.e., the LSB of scalar
+        ProjectivePoint::conditional_swap(&mut x0, &mut x1, Choice::from(prev_bit as u8));
+        // Don't leave the bit in the stack
+        #[cfg(feature = "zeroize")]
+        prev_bit.zeroize();
+
+        x0.as_affine()
+    }
+
     /// View this `MontgomeryPoint` as an array of bytes.
     pub const fn as_bytes(&self) -> &[u8; 32] {
         &self.0
@@ -357,55 +394,27 @@ define_mul_variants!(
 );
 
 /// Multiply this `MontgomeryPoint` by a `Scalar`.
-impl<'a, 'b> Mul<&'b Scalar> for &'a MontgomeryPoint {
+impl Mul<&Scalar> for &MontgomeryPoint {
     type Output = MontgomeryPoint;
 
-    /// Given `self` \\( = u\_0(P) \\), and a `Scalar` \\(n\\), return \\( u\_0(\[n\]P) \\).
-    fn mul(self, scalar: &'b Scalar) -> MontgomeryPoint {
-        // Algorithm 8 of Costello-Smith 2017
-        let affine_u = FieldElement::from_bytes(&self.0);
-        let mut x0 = ProjectivePoint::identity();
-        let mut x1 = ProjectivePoint {
-            U: affine_u,
-            W: FieldElement::ONE,
-        };
-
-        // NOTE: The below swap-double-add routine skips the first iteration, i.e., it assumes the
-        // MSB of `scalar` is 0. This is allowed, since it follows from Scalar invariant #1.
-
-        // Go through the bits from most to least significant, using a sliding window of 2
-        let mut bits = scalar.bits_le().rev();
-        let mut prev_bit = bits.next().unwrap();
-        for cur_bit in bits {
-            let choice: u8 = (prev_bit ^ cur_bit) as u8;
-
-            debug_assert!(choice == 0 || choice == 1);
-
-            ProjectivePoint::conditional_swap(&mut x0, &mut x1, choice.into());
-            differential_add_and_double(&mut x0, &mut x1, &affine_u);
-
-            prev_bit = cur_bit;
-        }
-        // The final value of prev_bit above is scalar.bits()[0], i.e., the LSB of scalar
-        ProjectivePoint::conditional_swap(&mut x0, &mut x1, Choice::from(prev_bit as u8));
-        // Don't leave the bit in the stack
-        #[cfg(feature = "zeroize")]
-        prev_bit.zeroize();
-
-        x0.as_affine()
+    /// Given `self` \\( = u\_0(P) \\), and a `Scalar` \\(n\\), return \\( u\_0(\[n\]P) \\)
+    fn mul(self, scalar: &Scalar) -> MontgomeryPoint {
+        // We multiply by the integer representation of the given Scalar. By scalar invariant #1,
+        // the MSB is 0, so we can skip it.
+        self.mul_bits_be(scalar.bits_le().rev().skip(1))
     }
 }
 
-impl<'b> MulAssign<&'b Scalar> for MontgomeryPoint {
-    fn mul_assign(&mut self, scalar: &'b Scalar) {
+impl MulAssign<&Scalar> for MontgomeryPoint {
+    fn mul_assign(&mut self, scalar: &Scalar) {
         *self = (self as &MontgomeryPoint) * scalar;
     }
 }
 
-impl<'a, 'b> Mul<&'b MontgomeryPoint> for &'a Scalar {
+impl Mul<&MontgomeryPoint> for &Scalar {
     type Output = MontgomeryPoint;
 
-    fn mul(self, point: &'b MontgomeryPoint) -> MontgomeryPoint {
+    fn mul(self, point: &MontgomeryPoint) -> MontgomeryPoint {
         point * self
     }
 }
@@ -422,7 +431,7 @@ mod test {
     #[cfg(feature = "alloc")]
     use alloc::vec::Vec;
 
-    use rand_core::RngCore;
+    use rand_core::{CryptoRng, RngCore};
 
     #[test]
     fn identity_in_different_coordinates() {
@@ -505,22 +514,99 @@ mod test {
         assert_eq!(u18, u18_unred);
     }
 
+    /// Returns a random point on the prime-order subgroup
+    fn rand_prime_order_point(mut rng: impl RngCore + CryptoRng) -> EdwardsPoint {
+        let s: Scalar = Scalar::random(&mut rng);
+        EdwardsPoint::mul_base(&s)
+    }
+
+    /// Given a bytestring that's little-endian at the byte level, return an iterator over all the
+    /// bits, in little-endian order.
+    fn bytestring_bits_le(x: &[u8]) -> impl DoubleEndedIterator<Item = bool> + Clone + '_ {
+        let bitlen = x.len() * 8;
+        (0..bitlen).map(|i| {
+            // As i runs from 0..256, the bottom 3 bits index the bit, while the upper bits index
+            // the byte. Since self.bytes is little-endian at the byte level, this iterator is
+            // little-endian on the bit level
+            ((x[i >> 3] >> (i & 7)) & 1u8) == 1
+        })
+    }
+
     #[test]
     fn montgomery_ladder_matches_edwards_scalarmult() {
         let mut csprng = rand_core::OsRng;
 
         for _ in 0..100 {
-            let s: Scalar = Scalar::random(&mut csprng);
-            let p_edwards = EdwardsPoint::mul_base(&s);
+            let p_edwards = rand_prime_order_point(&mut csprng);
             let p_montgomery: MontgomeryPoint = p_edwards.to_montgomery();
 
+            let s: Scalar = Scalar::random(&mut csprng);
             let expected = s * p_edwards;
             let result = s * p_montgomery;
 
             assert_eq!(result, expected.to_montgomery())
         }
     }
 
+    // Tests that, on the prime-order subgroup, MontgomeryPoint::mul_bits_be is the same as
+    // multiplying by the Scalar representation of the same bits
+    #[test]
+    fn montgomery_mul_bits_be() {
+        let mut csprng = rand_core::OsRng;
+
+        for _ in 0..100 {
+            // Make a random prime-order point P
+            let p_edwards = rand_prime_order_point(&mut csprng);
+            let p_montgomery: MontgomeryPoint = p_edwards.to_montgomery();
+
+            // Make a random integer b
+            let mut bigint = [0u8; 64];
+            csprng.fill_bytes(&mut bigint[..]);
+            let bigint_bits_be = bytestring_bits_le(&bigint).rev();
+
+            // Check that bP is the same whether calculated as scalar-times-edwards or
+            // integer-times-montgomery.
+            let expected = Scalar::from_bytes_mod_order_wide(&bigint) * p_edwards;
+            let result = p_montgomery.mul_bits_be(bigint_bits_be);
+            assert_eq!(result, expected.to_montgomery())
+        }
+    }
+
+    // Tests that MontgomeryPoint::mul_bits_be is consistent on any point, even ones that might be
+    // on the curve's twist. Specifically, this tests that b₁(b₂P) == b₂(b₁P) for random
+    // integers b₁, b₂ and random (curve or twist) point P.
+    #[test]
+    fn montgomery_mul_bits_be_twist() {
+        let mut csprng = rand_core::OsRng;
+
+        for _ in 0..100 {
+            // Make a random point P on the curve or its twist
+            let p_montgomery = {
+                let mut buf = [0u8; 32];
+                csprng.fill_bytes(&mut buf);
+                MontgomeryPoint(buf)
+            };
+
+            // Compute two big integers b₁ and b₂
+            let mut bigint1 = [0u8; 64];
+            let mut bigint2 = [0u8; 64];
+            csprng.fill_bytes(&mut bigint1[..]);
+            csprng.fill_bytes(&mut bigint2[..]);
+
+            // Compute b₁P and b₂P
+            let bigint1_bits_be = bytestring_bits_le(&bigint1).rev();
+            let bigint2_bits_be = bytestring_bits_le(&bigint2).rev();
+            let prod1 = p_montgomery.mul_bits_be(bigint1_bits_be.clone());
+            let prod2 = p_montgomery.mul_bits_be(bigint2_bits_be.clone());
+
+            // Check that b₁(b₂P) == b₂(b₁P)
+            assert_eq!(
+                prod1.mul_bits_be(bigint2_bits_be),
+                prod2.mul_bits_be(bigint1_bits_be)
+            );
+        }
+    }
+
     /// Check that mul_base_clamped and mul_clamped agree
     #[test]
     fn mul_base_clamped() {