Merge pull request #259 from dalek-cryptography/pippenger-bugfix

Pippenger bugfix

src/backend/serial/scalar_mul/pippenger.rs

@@ -88,14 +88,14 @@ impl VartimeMultiscalarMul for Pippenger {
         };
 
         let max_digit: usize = 1 << w;
-        let digits_count: usize = (256 + w - 1) / w; // == ceil(256/w)
+        let digits_count: usize = Scalar::to_radix_2w_size_hint(w);
         let buckets_count: usize = max_digit / 2; // digits are signed+centered hence 2^w/2, excluding 0-th bucket
 
         // Collect optimized scalars and points in buffers for repeated access
         // (scanning the whole set per digit position).
         let scalars = scalars
             .into_iter()
-            .map(|s| s.borrow().to_radix_2w(w).0);
+            .map(|s| s.borrow().to_radix_2w(w));
 
         let points = points
             .into_iter()
@@ -125,7 +125,8 @@ impl VartimeMultiscalarMul for Pippenger {
             // Note: if we add support for precomputed lookup tables,
             // we'll be adding/subtracting point premultiplied by `digits[i]` to buckets[0].
             for (digits, pt) in scalars_points.iter() {
-                let digit = digits[digit_index];
+                // Widen digit so that we don't run into edge cases when w=8.
+                let digit = digits[digit_index] as i16;
                 if digit > 0 {
                     let b = (digit - 1) as usize;
                     buckets[b] = (&buckets[b] + pt).to_extended();

src/backend/serial/scalar_mul/precomputed_straus.rs

@@ -85,7 +85,7 @@ impl VartimePrecomputedMultiscalarMul for VartimePrecomputedStraus {
         // nonzero NAF coefficient, but since we might have a lot of
         // them to search, it's not clear it's worthwhile to check.
         let mut S = ProjectivePoint::identity();
-        for j in (0..255).rev() {
+        for j in (0..256).rev() {
             let mut R: CompletedPoint = S.double();
 
             for i in 0..dp {

src/backend/serial/scalar_mul/straus.rs

@@ -179,7 +179,7 @@ impl VartimeMultiscalarMul for Straus {
 
         let mut r = ProjectivePoint::identity();
 
-        for i in (0..255).rev() {
+        for i in (0..256).rev() {
             let mut t: CompletedPoint = r.double();
 
             for (naf, lookup_table) in nafs.iter().zip(lookup_tables.iter()) {

src/backend/serial/scalar_mul/vartime_double_base.rs

@@ -23,7 +23,7 @@ pub fn mul(a: &Scalar, A: &EdwardsPoint, b: &Scalar) -> EdwardsPoint {
 
     // Find starting index
     let mut i: usize = 255;
-    for j in (0..255).rev() {
+    for j in (0..256).rev() {
         i = j;
         if a_naf[i] != 0 || b_naf[i] != 0 {
             break;

src/backend/vector/scalar_mul/pippenger.rs

@@ -45,14 +45,14 @@ impl VartimeMultiscalarMul for Pippenger {
         };
 
         let max_digit: usize = 1 << w;
-        let digits_count: usize = (256 + w - 1) / w; // == ceil(256/w)
+        let digits_count: usize = Scalar::to_radix_2w_size_hint(w);
         let buckets_count: usize = max_digit / 2; // digits are signed+centered hence 2^w/2, excluding 0-th bucket
 
         // Collect optimized scalars and points in a buffer for repeated access
         // (scanning the whole collection per each digit position).
         let scalars = scalars
             .into_iter()
-            .map(|s| s.borrow().to_radix_2w(w).0);
+            .map(|s| s.borrow().to_radix_2w(w));
 
         let points = points
             .into_iter()
@@ -82,7 +82,8 @@ impl VartimeMultiscalarMul for Pippenger {
             // Note: if we add support for precomputed lookup tables,
             // we'll be adding/subtractiong point premultiplied by `digits[i]` to buckets[0].
             for (digits, pt) in scalars_points.iter() {
-                let digit = digits[digit_index];
+                // Widen digit so that we don't run into edge cases when w=8.
+                let digit = digits[digit_index] as i16;
                 if digit > 0 {
                     let b = (digit - 1) as usize;
                     buckets[b] = &buckets[b] + pt;

src/backend/vector/scalar_mul/precomputed_straus.rs

@@ -84,7 +84,7 @@ impl VartimePrecomputedMultiscalarMul for VartimePrecomputedStraus {
         // nonzero NAF coefficient, but since we might have a lot of
         // them to search, it's not clear it's worthwhile to check.
         let mut R = ExtendedPoint::identity();
-        for j in (0..255).rev() {
+        for j in (0..256).rev() {
             R = R.double();
 
             for i in 0..dp {

src/backend/vector/scalar_mul/straus.rs

@@ -94,7 +94,7 @@ impl VartimeMultiscalarMul for Straus {
 
         let mut Q = ExtendedPoint::identity();
 
-        for i in (0..255).rev() {
+        for i in (0..256).rev() {
             Q = Q.double();
 
             for (naf, lookup_table) in nafs.iter().zip(lookup_tables.iter()) {

src/backend/vector/scalar_mul/vartime_double_base.rs

@@ -23,7 +23,7 @@ pub fn mul(a: &Scalar, A: &EdwardsPoint, b: &Scalar) -> EdwardsPoint {
 
     // Find starting index
     let mut i: usize = 255;
-    for j in (0..255).rev() {
+    for j in (0..256).rev() {
         i = j;
         if a_naf[i] != 0 || b_naf[i] != 0 {
             break;

src/edwards.rs

@@ -1255,6 +1255,73 @@ mod test {
         assert!(P1.compress().to_bytes() == P2.compress().to_bytes());
     }
 
+    // A single iteration of a consistency check for MSM.
+    fn multiscalar_consistency_iter(n: usize) {
+        use core::iter;
+        let mut rng = rand::thread_rng();
+
+        // Construct random coefficients x0, ..., x_{n-1},
+        // followed by some extra hardcoded ones.
+        let xs = (0..n)
+            .map(|_| Scalar::random(&mut rng))
+            // The largest scalar allowed by the type system, 2^255-1
+            .chain(iter::once(Scalar::from_bits([0xff; 32])))
+            .collect::<Vec<_>>();
+        let check = xs.iter()
+            .map(|xi| xi * xi)
+            .sum::<Scalar>();
+
+        // Construct points G_i = x_i * B
+        let Gs = xs.iter()
+            .map(|xi| xi * &constants::ED25519_BASEPOINT_TABLE)
+            .collect::<Vec<_>>();
+
+        // Compute H1 = <xs, Gs> (consttime)
+        let H1 = EdwardsPoint::multiscalar_mul(&xs, &Gs);
+        // Compute H2 = <xs, Gs> (vartime)
+        let H2 = EdwardsPoint::vartime_multiscalar_mul(&xs, &Gs);
+        // Compute H3 = <xs, Gs> = sum(xi^2) * B
+        let H3 = &check * &constants::ED25519_BASEPOINT_TABLE;
+
+        assert_eq!(H1, H3);
+        assert_eq!(H2, H3);
+    }
+
+    // Use different multiscalar sizes to hit different internal
+    // parameters.
+
+    #[test]
+    fn multiscalar_consistency_n_100() {
+        let iters = 50;
+        for _ in 0..iters {
+            multiscalar_consistency_iter(100);
+        }
+    }
+
+    #[test]
+    fn multiscalar_consistency_n_250() {
+        let iters = 50;
+        for _ in 0..iters {
+            multiscalar_consistency_iter(250);
+        }
+    }
+
+    #[test]
+    fn multiscalar_consistency_n_500() {
+        let iters = 50;
+        for _ in 0..iters {
+            multiscalar_consistency_iter(500);
+        }
+    }
+
+    #[test]
+    fn multiscalar_consistency_n_1000() {
+        let iters = 50;
+        for _ in 0..iters {
+            multiscalar_consistency_iter(1000);
+        }
+    }
+
     #[test]
     fn vartime_precomputed_vs_nonprecomputed_multiscalar() {
         let mut rng = rand::thread_rng();

src/scalar.rs

@@ -961,6 +961,24 @@ impl Scalar {
         output
     }
 
+    /// Returns a size hint indicating how many entries of the return
+    /// value of `to_radix_2w` are nonzero.
+    pub(crate) fn to_radix_2w_size_hint(w: usize) -> usize {
+        debug_assert!(w >= 6);
+        debug_assert!(w <= 8);
+
+        let digits_count = match w {
+            6 => (256 + w - 1)/w as usize,
+            7 => (256 + w - 1)/w as usize,
+            // See comment in to_radix_2w on handling the terminal carry.
+            8 => (256 + w - 1)/w + 1 as usize,
+            _ => panic!("invalid radix parameter"),
+        };
+
+        debug_assert!(digits_count <= 43);
+        digits_count
+    }
+
     /// Creates a representation of a Scalar in radix 64, 128 or 256 for use with the Pippenger algorithm.
     /// For lower radix, use `to_radix_16`, which is used by the Straus multi-scalar multiplication.
     /// Higher radixes are not supported to save cache space. Radix 256 is near-optimal even for very
@@ -979,13 +997,10 @@ impl Scalar {
     /// $$
     /// with \\(-2\^w/2 \leq a_i < 2\^w/2\\) for \\(0 \leq i < (n-1)\\) and \\(-2\^w/2 \leq a_{n-1} \leq 2\^w/2\\).
     ///
-    pub(crate) fn to_radix_2w(&self, w: usize) -> ([i8; 43], usize) {
+    pub(crate) fn to_radix_2w(&self, w: usize) -> [i8; 43] {
         debug_assert!(w >= 6);
         debug_assert!(w <= 8);
 
-        let digits_count = (256 + w - 1)/w as usize;
-        debug_assert!(digits_count <= 43);
-
         use byteorder::{ByteOrder, LittleEndian};
 
         // Scalar formatted as four `u64`s with carry bit packed into the highest bit.
@@ -997,6 +1012,7 @@ impl Scalar {
 
         let mut carry = 0u64;
         let mut digits = [0i8; 43];
+        let digits_count = (256 + w - 1)/w as usize;
         for i in 0..digits_count {
             // Construct a buffer of bits of the scalar, starting at `bit_offset`.
             let bit_offset = i*w;
@@ -1017,22 +1033,25 @@ impl Scalar {
             // Read the actual coefficient value from the window
             let coef = carry + (bit_buf & window_mask); // coef = [0, 2^r)
 
-             // Recenter coefficients from [0,2^r) to [-2^r/2, 2^r/2)
+             // Recenter coefficients from [0,2^w) to [-2^w/2, 2^w/2)
             carry = (coef + (radix/2) as u64) >> w;
             digits[i] = ((coef as i64) - (carry << w) as i64) as i8;
         }
 
-        // Apply the resulting carry to the last digit
-        // Since the highest bit of the 256-bit integer is 0,
-        // the last coefficient would always be in the lower half _inclusive_,
-        // so the carry in the end can be 1 iff the word equals 2^r/2.
-        // Since ±2^r/2 values are valid, to avoid adding an extra word,
-        // we allow the last word to touch the value 2^r/2.
-        // XXX: make sure tests cover this case, so the carry is non-zero and this line matters.
-        // Maybe it never happens to be non-zero for r=6/7/8?...
-        digits[digits_count-1] += (carry << w) as i8;
+        // When w < 8, we can fold the final carry onto the last digit d,
+        // because d < 2^w/2 so d + carry*2^w = d + 1*2^w < 2^(w+1) < 2^8.
+        //
+        // When w = 8, we can't fit carry*2^w into an i8.  This should
+        // not happen anyways, because the final carry will be 0 for
+        // reduced scalars, but the Scalar invariant allows 255-bit scalars.
+        // To handle this, we expand the size_hint by 1 when w=8,
+        // and accumulate the final carry onto another digit.
+        match w {
+            8 => digits[digits_count] += carry as i8,
+            _ => digits[digits_count-1] += (carry << w) as i8,
+        }
 
-        (digits, digits_count)
+        digits
     }
 
     /// Unpack this `Scalar` to an `UnpackedScalar` for faster arithmetic.
@@ -1512,32 +1531,43 @@ mod test {
         }
     }
 
+    fn test_pippenger_radix_iter(scalar: Scalar, w: usize) {
+        let digits_count = Scalar::to_radix_2w_size_hint(w);
+        let digits = scalar.to_radix_2w(w);
+
+        let radix = Scalar::from((1<<w) as u64);
+        let mut term = Scalar::one();
+        let mut recovered_scalar = Scalar::zero();
+        for digit in &digits[0..digits_count] {
+            let digit = *digit;
+            if digit != 0 {
+                let sdigit = if digit < 0 {
+                    -Scalar::from((-(digit as i64)) as u64)
+                } else {
+                    Scalar::from(digit as u64)
+                };
+                recovered_scalar += term * sdigit;
+            }
+            term *= radix;
+        }
+        // When the input is unreduced, we may only recover the scalar mod l.
+        assert_eq!(recovered_scalar, scalar.reduce());
+    }
+
     #[test]
     fn test_pippenger_radix() {
         use core::iter;
         // For each valid radix it tests that 1000 random-ish scalars can be restored
         // from the produced representation precisely.
-        for w in 6..9 {
-            for scalar in (2..100).map(|s| Scalar::from(s as u64).invert() ).chain(iter::once(-Scalar::one())) {
-                let (digits, digits_count) = scalar.to_radix_2w(w);
+        let cases = (2..100)
+            .map(|s| Scalar::from(s as u64).invert())
+            // The largest unreduced scalar, s = 2^255-1
+            .chain(iter::once(Scalar::from_bits([0xff; 32])));
 
-                let radix = Scalar::from((1<<w) as u64);
-                let mut term = Scalar::one();
-                let mut recovered_scalar = Scalar::zero();
-                for digit in &digits[0..digits_count] {
-                    let digit = *digit;
-                    if digit != 0 {
-                        let sdigit = if digit < 0 {
-                            -Scalar::from((-(digit as i64)) as u64)
-                        } else {
-                            Scalar::from(digit as u64)
-                        };
-                        recovered_scalar += term * sdigit;
-                    }
-                    term *= radix;
-                }
-                assert_eq!(recovered_scalar, scalar);
-            }
+        for scalar in cases {
+            test_pippenger_radix_iter(scalar, 6);
+            test_pippenger_radix_iter(scalar, 7);
+            test_pippenger_radix_iter(scalar, 8);
         }
     }
 }