Merge pull request #249 from oleganza/oleg/pippenger2

Pippenger multiscalar multiplication algorithm

src/backend/serial/scalar_mul/mod.rs

@@ -26,3 +26,6 @@ pub mod straus;
 
 #[cfg(feature = "alloc")]
 pub mod precomputed_straus;
+
+#[cfg(feature = "alloc")]
+pub mod pippenger;

src/backend/serial/scalar_mul/pippenger.rs

@@ -0,0 +1,205 @@
+// -*- mode: rust; -*-
+//
+// This file is part of curve25519-dalek.
+// Copyright (c) 2019 Oleg Andreev
+// See LICENSE for licensing information.
+//
+// Authors:
+// - Oleg Andreev <oleganza@gmail.com>
+
+//! Implementation of a variant of Pippenger's algorithm.
+
+#![allow(non_snake_case)]
+
+use core::borrow::Borrow;
+
+use edwards::EdwardsPoint;
+use scalar::Scalar;
+use traits::VartimeMultiscalarMul;
+
+#[allow(unused_imports)]
+use prelude::*;
+
+/// Implements a version of Pippenger's algorithm.
+///
+/// The algorithm works as follows:
+///
+/// Let `n` be a number of point-scalar pairs.
+/// Let `w` be a window of bits (6..8, chosen based on `n`, see cost factor).
+///
+/// 1. Prepare `2^(w-1) - 1` buckets with indices `[1..2^(w-1))` initialized with identity points.
+///    Bucket 0 is not needed as it would contain points multiplied by 0.
+/// 2. Convert scalars to a radix-`2^w` representation with signed digits in `[-2^w/2, 2^w/2]`.
+///    Note: only the last digit may equal `2^w/2`.
+/// 3. Starting with the last window, for each point `i=[0..n)` add it to a a bucket indexed by
+///    the point's scalar's value in the window.
+/// 4. Once all points in a window are sorted into buckets, add buckets by multiplying each
+///    by their index. Efficient way of doing it is to start with the last bucket and compute two sums:
+///    intermediate sum from the last to the first, and the full sum made of all intermediate sums.
+/// 5. Shift the resulting sum of buckets by `w` bits by using `w` doublings.
+/// 6. Add to the return value.
+/// 7. Repeat the loop.
+///
+/// Approximate cost w/o wNAF optimizations (A = addition, D = doubling):
+///
+/// ```ascii
+/// cost = (n*A + 2*(2^w/2)*A + w*D + A)*256/w
+///          |          |       |     |   |
+///          |          |       |     |   looping over 256/w windows
+///          |          |       |     adding to the result
+///    sorting points   |       shifting the sum by w bits (to the next window, starting from last window)
+///    one by one       |
+///    into buckets     adding/subtracting all buckets
+///                     multiplied by their indexes
+///                     using a sum of intermediate sums
+/// ```
+///
+/// For large `n`, dominant factor is (n*256/w) additions.
+/// However, if `w` is too big and `n` is not too big, then `(2^w/2)*A` could dominate.
+/// Therefore, the optimal choice of `w` grows slowly as `n` grows.
+///
+/// This algorithm is adapted from section 4 of https://eprint.iacr.org/2012/549.pdf.
+pub struct Pippenger;
+
+#[cfg(any(feature = "alloc", feature = "std"))]
+impl VartimeMultiscalarMul for Pippenger {
+    type Point = EdwardsPoint;
+
+    fn optional_multiscalar_mul<I, J>(scalars: I, points: J) -> Option<EdwardsPoint>
+    where
+        I: IntoIterator,
+        I::Item: Borrow<Scalar>,
+        J: IntoIterator<Item = Option<EdwardsPoint>>,
+    {
+        use traits::Identity;
+
+        let mut scalars = scalars.into_iter();
+        let size = scalars.by_ref().size_hint().0;
+
+        // Digit width in bits. As digit width grows,
+        // number of point additions goes down, but amount of
+        // buckets and bucket additions grows exponentially.
+        let w = if size < 500 {
+            6
+        } else if size < 800 {
+            7
+        } else {
+            8
+        };
+
+        let max_digit: usize = 1 << w;
+        let digits_count: usize = (256 + w - 1) / w; // == ceil(256/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);
+
+        let points = points
+            .into_iter()
+            .map(|p| p.map(|P| P.to_projective_niels()));
+
+        let scalars_points = scalars.zip(points).map(|(s,maybe_p)| maybe_p.map(|p| (s,p) ) )
+            .collect::<Option<Vec<_>>>();
+        let scalars_points = match scalars_points {
+            Some(sp) => sp,
+            None => return None,
+        };
+
+        // Prepare 2^w/2 buckets.
+        // buckets[i] corresponds to a multiplication factor (i+1).
+        let mut buckets: Vec<_> = (0..buckets_count)
+            .map(|_| EdwardsPoint::identity())
+            .collect();
+
+        let mut columns = (0..digits_count).rev().map(|digit_index| {
+            // Clear the buckets when processing another digit.
+            for i in 0..buckets_count {
+                buckets[i] = EdwardsPoint::identity();
+            }
+
+            // Iterate over pairs of (point, scalar)
+            // and add/sub the point to the corresponding bucket.
+            // 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];
+                if digit > 0 {
+                    let b = (digit - 1) as usize;
+                    buckets[b] = (&buckets[b] + pt).to_extended();
+                } else if digit < 0 {
+                    let b = (-digit - 1) as usize;
+                    buckets[b] = (&buckets[b] - pt).to_extended();
+                }
+            }
+
+            // Add the buckets applying the multiplication factor to each bucket.
+            // The most efficient way to do that is to have a single sum with two running sums:
+            // an intermediate sum from last bucket to the first, and a sum of intermediate sums.
+            //
+            // For example, to add buckets 1*A, 2*B, 3*C we need to add these points:
+            //   C
+            //   C B
+            //   C B A   Sum = C + (C+B) + (C+B+A)
+            let mut buckets_intermediate_sum = buckets[buckets_count - 1];
+            let mut buckets_sum = buckets[buckets_count - 1];
+            for i in (0..(buckets_count - 1)).rev() {
+                buckets_intermediate_sum += buckets[i];
+                buckets_sum += buckets_intermediate_sum;
+            }
+
+            buckets_sum
+        });
+
+        // Take the high column as an initial value to avoid wasting time doubling the identity element in `fold()`.
+        // `unwrap()` always succeeds because we know we have more than zero digits.
+        let hi_column = columns.next().unwrap();
+
+        Some(
+            columns
+                .fold(hi_column, |total, p| total.mul_by_pow_2(w as u32) + p)
+                .into(),
+        )
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use constants;
+    use scalar::Scalar;
+
+    #[test]
+    fn test_vartime_pippenger() {
+        // Reuse points across different tests
+        let mut n = 512;
+        let x = Scalar::from(2128506u64).invert();
+        let y = Scalar::from(4443282u64).invert();
+        let points: Vec<_> = (0..n)
+            .map(|i| constants::ED25519_BASEPOINT_POINT * Scalar::from(1 + i as u64))
+            .collect();
+        let scalars: Vec<_> = (0..n)
+            .map(|i| x + (Scalar::from(i as u64) * y)) // fast way to make ~random but deterministic scalars
+            .collect();
+
+        let premultiplied: Vec<EdwardsPoint> = scalars
+            .iter()
+            .zip(points.iter())
+            .map(|(sc, pt)| sc * pt)
+            .collect();
+
+        while n > 0 {
+            let scalars = &scalars[0..n].to_vec();
+            let points = &points[0..n].to_vec();
+            let control: EdwardsPoint = premultiplied[0..n].iter().sum();
+
+            let subject = Pippenger::vartime_multiscalar_mul(scalars.clone(), points.clone());
+
+            assert_eq!(subject.compress(), control.compress());
+
+            n = n / 2;
+        }
+    }
+}

src/backend/vector/scalar_mul/mod.rs

@@ -17,3 +17,6 @@ pub mod straus;
 
 #[cfg(feature = "alloc")]
 pub mod precomputed_straus;
+
+#[cfg(feature = "alloc")]
+pub mod pippenger;

src/backend/vector/scalar_mul/pippenger.rs

@@ -0,0 +1,165 @@
+// -*- mode: rust; -*-
+//
+// This file is part of curve25519-dalek.
+// Copyright (c) 2019 Oleg Andreev
+// See LICENSE for licensing information.
+//
+// Authors:
+// - Oleg Andreev <oleganza@gmail.com>
+
+#![allow(non_snake_case)]
+
+use core::borrow::Borrow;
+
+use backend::vector::{CachedPoint, ExtendedPoint};
+use edwards::EdwardsPoint;
+use scalar::Scalar;
+use traits::{Identity, VartimeMultiscalarMul};
+
+#[allow(unused_imports)]
+use prelude::*;
+
+/// Implements a version of Pippenger's algorithm.
+///
+/// See the documentation in the serial `scalar_mul::pippenger` module for details.
+pub struct Pippenger;
+
+#[cfg(any(feature = "alloc", feature = "std"))]
+impl VartimeMultiscalarMul for Pippenger {
+    type Point = EdwardsPoint;
+
+    fn optional_multiscalar_mul<I, J>(scalars: I, points: J) -> Option<EdwardsPoint>
+    where
+        I: IntoIterator,
+        I::Item: Borrow<Scalar>,
+        J: IntoIterator<Item = Option<EdwardsPoint>>,
+    {
+        let mut scalars = scalars.into_iter();
+        let size = scalars.by_ref().size_hint().0;
+        let w = if size < 500 {
+            6
+        } else if size < 800 {
+            7
+        } else {
+            8
+        };
+
+        let max_digit: usize = 1 << w;
+        let digits_count: usize = (256 + w - 1) / w; // == ceil(256/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);
+
+        let points = points
+            .into_iter()
+            .map(|p| p.map(|P| CachedPoint::from(ExtendedPoint::from(P))));
+
+        let scalars_points = scalars.zip(points).map(|(s,maybe_p)| maybe_p.map(|p| (s,p) ) )
+            .collect::<Option<Vec<_>>>();
+        let scalars_points = match scalars_points {
+            Some(sp) => sp,
+            None => return None,
+        };
+
+        // Prepare 2^w/2 buckets.
+        // buckets[i] corresponds to a multiplication factor (i+1).
+        let mut buckets: Vec<ExtendedPoint> = (0..buckets_count)
+            .map(|_| ExtendedPoint::identity())
+            .collect();
+
+        let mut columns = (0..digits_count).rev().map(|digit_index| {
+            // Clear the buckets when processing another digit.
+            for i in 0..buckets_count {
+                buckets[i] = ExtendedPoint::identity();
+            }
+
+            // Iterate over pairs of (point, scalar)
+            // and add/sub the point to the corresponding bucket.
+            // 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];
+                if digit > 0 {
+                    let b = (digit - 1) as usize;
+                    buckets[b] = &buckets[b] + pt;
+                } else if digit < 0 {
+                    let b = (-digit - 1) as usize;
+                    buckets[b] = &buckets[b] - pt;
+                }
+            }
+
+            // Add the buckets applying the multiplication factor to each bucket.
+            // The most efficient way to do that is to have a single sum with two running sums:
+            // an intermediate sum from last bucket to the first, and a sum of intermediate sums.
+            //
+            // For example, to add buckets 1*A, 2*B, 3*C we need to add these points:
+            //   C
+            //   C B
+            //   C B A   Sum = C + (C+B) + (C+B+A)
+            let mut buckets_intermediate_sum = buckets[buckets_count - 1];
+            let mut buckets_sum = buckets[buckets_count - 1];
+            for i in (0..(buckets_count - 1)).rev() {
+                buckets_intermediate_sum =
+                    &buckets_intermediate_sum + &CachedPoint::from(buckets[i]);
+                buckets_sum = &buckets_sum + &CachedPoint::from(buckets_intermediate_sum);
+            }
+
+            buckets_sum
+        });
+
+        // Take the high column as an initial value to avoid wasting time doubling the identity element in `fold()`.
+        // `unwrap()` always succeeds because we know we have more than zero digits.
+        let hi_column = columns.next().unwrap();
+
+        Some(
+            columns
+                .fold(hi_column, |total, p| {
+                    &total.mul_by_pow_2(w as u32) + &CachedPoint::from(p)
+                })
+                .into(),
+        )
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use constants;
+    use scalar::Scalar;
+
+    #[test]
+    fn test_vartime_pippenger() {
+        // Reuse points across different tests
+        let mut n = 512;
+        let x = Scalar::from(2128506u64).invert();
+        let y = Scalar::from(4443282u64).invert();
+        let points: Vec<_> = (0..n)
+            .map(|i| constants::ED25519_BASEPOINT_POINT * Scalar::from(1 + i as u64))
+            .collect();
+        let scalars: Vec<_> = (0..n)
+            .map(|i| x + (Scalar::from(i as u64) * y)) // fast way to make ~random but deterministic scalars
+            .collect();
+
+        let premultiplied: Vec<EdwardsPoint> = scalars
+            .iter()
+            .zip(points.iter())
+            .map(|(sc, pt)| sc * pt)
+            .collect();
+
+        while n > 0 {
+            let scalars = &scalars[0..n].to_vec();
+            let points = &points[0..n].to_vec();
+            let control: EdwardsPoint = premultiplied[0..n].iter().sum();
+
+            let subject = Pippenger::vartime_multiscalar_mul(scalars.clone(), points.clone());
+
+            assert_eq!(subject.compress(), control.compress());
+
+            n = n / 2;
+        }
+    }
+}

src/edwards.rs

@@ -676,11 +676,15 @@ impl VartimeMultiscalarMul for EdwardsPoint {
         assert_eq!(s_hi, Some(s_lo));
         assert_eq!(p_hi, Some(p_lo));
 
-        // Now we know there's a single size.  When we do
+        // Now we know there's a single size.
-        // size-dependent algorithm dispatch, use this as the hint.
+        // Use this as the hint to decide which algorithm to use.
-        let _size = s_lo;
+        let size = s_lo;
 
-        scalar_mul::straus::Straus::optional_multiscalar_mul(scalars, points)
+        if size < 190 {
+            scalar_mul::straus::Straus::optional_multiscalar_mul(scalars, points)
+        } else {
+            scalar_mul::pippenger::Pippenger::optional_multiscalar_mul(scalars, points)
+        }
     }
 }
 

src/scalar.rs

@@ -961,6 +961,80 @@ impl Scalar {
         output
     }
 
+    /// 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
+    /// large inputs.
+    ///
+    /// Radix below 64 or above 256 is prohibited.
+    /// This method returns digits in a fixed-sized array, excess digits are zeroes.
+    /// The second returned value is the number of digits.
+    ///
+    /// ## Scalar representation
+    ///
+    /// Radix \\(2\^w\\), with \\(n = ceil(256/w)\\) coefficients in \\([-(2\^w)/2,(2\^w)/2)\\),
+    /// i.e., scalar is represented using digits \\(a\_i\\) such that
+    /// $$
+    ///    a = a\_0 + a\_1 2\^1w + \cdots + a_{n-1} 2\^{w*(n-1)},
+    /// $$
+    /// 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) {
+        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.
+        let mut scalar64x4 = [0u64; 4];
+        LittleEndian::read_u64_into(&self.bytes, &mut scalar64x4[0..4]);
+
+        let radix: u64 = 1 << w;
+        let window_mask: u64 = radix - 1;
+
+        let mut carry = 0u64;
+        let mut digits = [0i8; 43];
+        for i in 0..digits_count {
+            // Construct a buffer of bits of the scalar, starting at `bit_offset`.
+            let bit_offset = i*w;
+            let u64_idx = bit_offset / 64;
+            let bit_idx = bit_offset % 64;
+
+            // Read the bits from the scalar
+            let bit_buf: u64;
+            if bit_idx < 64 - w  || u64_idx == 3 {
+                // This window's bits are contained in a single u64,
+                // or it's the last u64 anyway.
+                bit_buf = scalar64x4[u64_idx] >> bit_idx;
+            } else {
+                // Combine the current u64's bits with the bits from the next u64
+                bit_buf = (scalar64x4[u64_idx] >> bit_idx) | (scalar64x4[1+u64_idx] << (64 - bit_idx));
+            }
+
+            // 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)
+            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;
+
+        (digits, digits_count)
+    }
+
     /// Unpack this `Scalar` to an `UnpackedScalar` for faster arithmetic.
     pub(crate) fn unpack(&self) -> UnpackedScalar {
         UnpackedScalar::from_bytes(&self.bytes)
@@ -1437,4 +1511,33 @@ mod test {
             assert_eq!(a * b, Scalar::one());
         }
     }
+
+    #[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 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);
+            }
+        }
+    }
 }